Difference between revisions of "Task 2: Alignments"

Latest revision as of 13:38, 1 September 2013

Sequence alignments are a good start to analyse a new protein sequence. Alignment search tools such as Blast, Psiblast or HHblits can be used to find related protein sequences.

Lab journal task 2

Sequence Searches

The HFE protein sequence was used to conduct sequence searches with the three different search tools Blast, Psiblast and HHblits. CATH and COPS, two protein structure classification databases, were used to evaluate and compare the results of the different tools.

CATH is a database that hierarchically classifies all domain structures from PDB. This hierarchy consists of four major levels: 1. Class 2. Architecture 3. Topology 4. Homologous superfamily

We compared the first level of the CATH hierarchy (Class) to check whether the hits have the same fold class as the query. The class of a domain is determined by its secondary structure composition (mostly alpha, mostly beta, mixed alpha/beta or few secondary structures). Each domain of a protein is assigned a CATH class, this means, that one query protein can consist of several fold classes, one for each domain. The HFE protein consists of two domains that belong to the CATH fold class 2 and 3.

Another database for protein structure clustering is COPS. It consists of different structural groups that contain all proteins with a structural similarity above a defined threshold. For example, all proteins in the L30 group have a structural similarity of at least 30% percent. This means that if two proteins are not in the same L30 group, they have less than 30% structural similarity. The L30 groups of the hits were determined to check if they are in the same group as the HFE protein or not.

Blast

The Blast search in the big_80 database yielded 1504 hits.

**Figure 1:** E-value and identity distribution of the Blast search hits.
a) E-value distribution of the results of the Blast search against big_80 using standard parameters. The x-axis shows the base 10 logarithm of the e-value and the y-axis the frequency. The better the e-value the smaller the log(e-value).	b) Identity distribution of the results of the Blast search against big_80 using standard parameters. The x-axis shows the relative identity between the query and the hits.

</figure>

<xr id="blast distribution"/> shows the e-value and identity distributions of all unfiltered Blast search hits. There are only a few hits with a very low e-value of 1e-160, but most significant hits have an e-value between 1e-40 and 1e-50 (left small peak in a)). Most results have a worse e-value above 1e-20. The histogram of the sequence identity b) shows that nearly all hits have an identity between 20% to 40% to the query.

The Blast results only contain 13 PDB IDs. Nevertheless, the number of shared CATH classes between the reference sequence HFE and the hits were computed as described in the lab journal.

<css> table.colBasic2 { margin-left: auto; margin-right: auto; border: 2px solid black; border-collapse:collapse; width: 40%; } .colBasic2 th,td { padding: 3px; border: 2px solid black; } .colBasic2 td { text-align:left; } .colBasic2 tr th { background-color:#efefef; color: black;} .colBasic2 tr:first-child th { background-color:#adceff; color:black;} </css>

**Table 1:** Histogram of the number of shared CATH classes between HFE and the Blast hits.
number of same fold classes	frequency
0	4
1	6
2	3

</figtable>

<xr id="blast cath"/> shows, that only 3 of the 13 PDB hits share all domain fold classes with HFE, 6 of the hits have one domain fold class in commom with the reference and 4 even none. Besides, only 6 of the 13 pdb hits are in the same COPS L30 group as the HFE protein. Those numbers indicate, that some of the Blast hits are not really related to HFE and thus, false positives.

Psiblast

Four different Psiblast runs against big_80 were executed using the following four parameter combinations:

2 iterations and default e-value cutoff 0.002 (j=2, h=2)
2 iterations and e-value cutoff 10E-10 (j=2, h=10E-10)
10 iterations and default e-value cutoff 0.002 (j=10, h=2)
10 iterations and e-value cutoff 10E-10 (j=10, h=10E-10)

The e-value and sequence identity distributions of the hits from the four runs were compared to identify the best parameter combination.

**Figure 2:** E-value and percent sequence identity distribution of the results from the four Psiblast runs against big_80.
a) Boxplot of the log(e-value) of the four different Psiblast searches against big_80: 2 iterations with default e-value cutoff (0.002) and 10E-10 and 10 iterations with default e-value cutoff (0.002) and 10E-10. The thick black line inside the box is the median and the upper and lower end of the box the 25% and 75% quantils. The single points above the whiskers are outliers.	b) Boxplot of the identity of the results of the four different Psiblast searches against big_80: 2 iterations with default e-value cutoff (0.002) and 10E-10 and 10 iterations with default e-value cutoff (0.002) and 10E-10. The thick black line inside the box is the median and the upper and lower end of the box the 25% and 75% quantils. The single points above and below the whiskers are outliers.

</figure>

**Table 2:** Comparison of the results of the four Psiblast searches against big_80.
		Psiblast parameters
		j=2, h=0.002	j=2, h=10E-10	j=10, h=0.002	j=10, h=10E-10
# hits		2564	3163	2601	3343
# PDB hits		66	72	56	66
CATH classes in common with HFE	0	21	21	20	21
	1	42	48	33	42
	2	3	3	3	3
# hits within same L30 group		6	6	6	6

</figtable>

<xr id="psiblast summary"/> summarises the results of the different Psiblast runs. The run with default parameter settings (j=2 and h=0.002) found the least hits. Nevertheless, the boxplot of the e-value distribution (<xr id="psiblast boxplot"/> a)) shows, that the default parameters can be considered the best: the median of the log(e-value) is smaller than the median of the other 3 psiblast runs. This indicates that the percentage of hits with a low and thus significant e-value is higher in the default psiblast run than in the other three. Besides, the default parameters also led to a higher percent sequence identity as indicated by the median of the cyan box in <xr id="psiblast boxplot"/> b).

The number of PDB hits varies among the different runs, but all found 3 proteins with the same two CATH classes as the query. There are 1095 proteins in the same COPS L30 group as the query HFE. Since the big_80 database does not contain all proteins from the PDB, the quality of the program cannot be assessed using the COPS L30 group of HFE. Nevertheless, from all found PDB hits, only 6 proteins belonged to the same L30 group.

PDB search

To get more PDB ids, we also did 4 Psiblast searches against the PDB seqres database using the profiles created with big_80. Those four result lists were compared as well.

**Figure 3:** E-value and sequence identity distributions of the hits from the four Psiblast runs against PDB_seqres.
b) Boxplot of the log(e-value) of the four different Psiblast searches against PDB_seqres using the checkpoint files from the earlier runs. (2 iterations with default e-value cutoff (0.002) and 10E-10 and 10 iterations with default e-value cutoff(0.002) and 10E-10). The thick black line inside the box is the median and the upper and lower end of the box the 25% and 75% quantils. The single points below the whiskers are outliers.	Boxplot of the identity of the results of the four different Psiblast searches against pdb_seqres using the checkpoint files from the runs before. (2 iterations with default e-value cutoff (0.002) and 10E-10 and 10 iterations with default e-value cutoff(0.002) and 10E-10). The thick black line inside the box is the median and the upper and lower end of the box the 25% and 75% quantils. The single points above and below the whiskers are outliers.

</figure>

**Table 3:** Comparison of the results of the four Psiblast searches against PDB_seqres.
		Psiblast parameters
		j=2, h=0.002	j=2, h=10E-10	j=10, h=0.002	j=10, h=10E-10
# hits		4307	4814	3874	4058
CATH classes in common with HFE	0	561	586	529	539
	1	2798	3280	2397	2571
	2	948	948	948	948
# hits within same L30 group		1047	1053	1045	1052

</figtable>

The e-value and sequence identity distributions in <xr id="psiblast pdb"/> do not show a significant difference in the results, because the medians of all four e-value distributions are nearly the same. The last run with the parameter combination of j=10 and h=10E-10 has a higher percentage of low e-value hits, because the 25% quantil is below those of the other distributions. But the run with j=2 and h=10E-10 lead to more hits with a higher sequence identity, because the median slightly lies above the other medians. Apart from the query sequence, there were no PDB hits with a sequence identity above 40%.

Evaluating the CATH and SCOP annotations of the found hits revealed that all four searches found 948 hits with the same 2 CATH classes as reference. Besides, all runs could find nearly all 1095 proteins that are in the same COPS L30 group as HFE. Therefore, the decision which of the parameter settings leads to the best results cannot be based on the PDB hits, since the differences are minimal. Nevertheless, j=10, h=0.002 had the least hits, but the highest percentage of hits in the same L30 group as the query. It is thus best suited to find a high number of significant hits.

**Figure 4:** Venn diagram of the overlap between all hits from the four different searches. J1H1: 10 iterations and e-value cutoff 10E-10, J1H2: 10 iterations and e-value cutoff 2E-3, J2H1: 2 iterations and e-value cutoff 10E-10, J2H2: 2 iterations and e-value cutoff 2E-3.

</figure>

We also computed the hits overlap of the four different Psiblast searches in big_80, to see how comparable the found hits are. The overlap is illustrated in the Venn diagram in <xr id="psiblast venn"/>.

Additionally, we computed the overlap for the 100 best hits (based on e-value) of each of the runs (data not shown). They were exactly the same in all four runs. Thus, the parameter settings only matter if a high number of hits is needed, but not if onyl the best hits of a search are used in further analyses.

HHblits

We executed HHblits with default parameters and searched the uniprot20_02Sep11 database. HHblits is an alignment search tool that searches in a database of clustered sequences. A profile is generated from the query that is compared with the profile of each cluster. For each cluster, not only the representative, but also all other cluster members are reported in the results list. Thus, HHblits finds much more sequences than Blast and Psiblast. We found 2861 clusters and 37908 sequences in total. We used all sequences for the analysis, because one of HHblits' most important features is the high number of good results.

**Figure 5:** E-value and identity distributions of the HHblits search hits.
a) E-value distribution of the results of the HHblits search with default parameters using the uniprot20_02Sep11 database. The x-axis shows the base 10 logarithm of the e-value and the y-axis the frequency. A smaller log10(e-value) corresponds to a better results.	b) Identity distribution of the results of the HHblits search in the uniprot20_02Sep11 database using default parameters. The x-axis shows the relative identity between the query and the hits.

</figure>

<xr id ="hhblits distribution"/> a) shows the distribution of the log10 of the e-values of the HHblits hits. It is remarkable that HHblits finds so many hits with a very low e-value of 1E-130. One reason is, that the best hit belongs to a very large cluster. Nevertheless, in comparison to the Blast e-value distribution ( <xr id="blast distribution"/> a)), the percentage of hits with a high e-value is much lower in the HHblits results.

There are several hits with a very low sequence identity to the query (see <xr id="hhblits distribution"/> b)), but the distribution contains two peaks around 20% and 40% identity that correspond to distant and more closely related hits.

Comparison

**Figure 6:** E-value and identity distribution of the Blast, Psiblast and HHblits search hits.
a) E-value distribution of the Blast, Psiblast and HHblits hits with default parameter searches.	b) Identity distribution of the Blast, Psiblast and HHblits hits default parameters searches

</figure>

We used the results of the Psiblast run with 2 iterations and an e-value cutoff of 2e-3, because we also used default parameters for Blast and HHblits and because we also received the best results with this parameters. We compared the e-value and identity distributions. <xr id ="comparison"/> a) clearly shows that Blast has the lowest e-values, but the percentage of low e-value hits is higher in the HHblits results. The median sequence identity has the highest value among the Blast hits.

**Figure 7:** Venn diagram of the overlap between all hits from the Blast, Psiblast and HHblits runs with default parameters.

</figure>

The overlap between all hits and the top 100 hits are shown in <xr id ="comparison venn"/>.

**Table 4:** Rank of the top 6 Psiblast hits in the Blast and HHblits results.
search tool	rank of hit
Psiblast	1	2	3	4	5	6
Blast	1	3	2	52	109	81
HHblits	2	not found	not found	not found	not found	not found

</figtable>

We tried to identify the 6 best Psiblast hits with the lowest e-values in the results of the Blast and HHblits searches, see <xr id ="best hits"/>- The top 6 Psiblast hits are all present in the Blast results and the first 3 hits are the same. This confirms that Blast and Psiblast are very similar. HHblits, on the contrary, has definitely not the same results as Blast and Psiblast. From the top 6 Psiblast hits, only the first one is present in the second best cluster. This proves, that HHblits finds more distant related hits, which is also indicated by the low median sequence identity of the HHblits hits.

The best 5 hits from all three search tools belong to the MHC class I family and are thus correct hits.

Multiple sequence alignments

**Table 5:** List of all sequences that were used in the multiple sequence alignments.
< 30% sequence identity		> 60% sequence identity		whole range of sequence identity
AC Number	Sequence Identity	AC Number	Sequence Identity	AC Number	Sequence Identity	AC Number	Sequence Identity
F6JYA9	0.29	B4DDZ1	0.83	B4DDZ1	0.83	F6WCX4	0.34
D9J389	0.28	Q5EEZ1	0.76	Q5EEZ1	0.76	F6JYA9	0.29
D5MSB3	0.27	G3THV5	0.75	G1MBW1	0.73	D9J389	0.28
B3FRK2	0.20	G1MBW1	0.73	H0VAR7	0.72	D5MSB3	0.27
3ov6_A	0.22	H0VAR7	0.72	F1PX48	0.71	2zok_E	0.24
1p7k_L	0.21	F1PX48	0.71	G1T7D7	0.70	H0Y1D0	0.20
H0Y1D0	0.20	G1T7D7	0.70	G5BQE5	0.67	Q8SNJ4	0.22
B3FRK3	0.18	O35799	0.68	Q95IT9	0.38	B3FRK3	0.18
Q8HWL2	0.17	G5BQE5	0.67	2qrt_A	0.38	Q8HWL2	0.17
				Q8HX83	0.36

In order to assess the difference of multiple sequence alignments between close and related homologs, three different groups of sequences were selected. One containing only sequences with a sequence identity above 60% to human hfe, one with sequences below 30% identity and one with sequences covering the whole range of sequence identity. The selected sequences are listed in Table 5.
In order to ensure that also sequences with known structures are included in the alignments, the sequences from the following PDB structures were selected:

1p7k_L
3ov6_A
2qrt_A
2zok_E

ClustalW

**Figure 8:** ClustalW alignments for all different groups of sequence identity are shown.
Below 30% identity.	Above 60% identity.	Whole range of identity.

</figtable>

**Figure 9:** The multiple alignments are displayed with the secondary structure annotation of Q30201 from Uniprot shown in red (helix), blue (sheet) and green (loop).
ClustalW alignment of sequences from the below 30% identity group.	ClustalW alignment of sequences from the above 60% identity group.	ClustalW alignment of sequences from the whole range identity group.

</figtable>

MAFFT

**Figure 10:** MAFFT alignments for all different groups of sequence identity are shown.
Below 30% identity.	Above 60% identity.	Whole range of identity.

</figtable>

T-Coffee

**Figure 11:** T-Coffee alignments for all different groups of sequence identity are shown.
Below 30% identity.	Above 60% identity.	Whole range of identity.

</figtable>

Expresso

**Figure 12:** Expresso alignments for all different groups of sequence identity are shown.
Below 30% identity.	Above 60% identity.	Whole range of identity.

</figtable>

Discussion

In the above 60% sequence identity group, the residues and the gaps are well conserved, especially in the first third of the sequence. In the other two thirds, the sequence conservation drops slightly, but is still at a generally high level. In the below 30% identity group, the number of gaps inside the sequences is not much higher than in the above 60% group, but the overall residue conservation is significantly lower. Nevertheless, there are some very well conserved residues in this group that might be functionally important. In the mixed sequence identity group, the gaps are not as well conserved as in the other two groups and the conservation is even lower than in the below 30% identity group. But this effect is probably strengthened by the higher amount of sequences present. Another peculiarity is that the inclusion of 3D information with Expresso (figure 12), yields an alignment with a lot more gaps, but more conserved columns than the other programs. Although this could be a desirable property in an alignment, this behavior can be reproduced with the other programs by modifying the gap penalty parameters.

The different alignment programs yield comparable results, i.e. none of them is considerably better or worse than the other two. Nevertheless, some differences can be observed. The first notable one is that MAFFT yields a lot less conserved columns for the low sequence identity group than the other programs. Also, the positioning of consecutive gap columns varies between the programs.

Since there are only one to two bigger gaps in the alignments and the secondary structure segments are often not very long, it is hard to tell if the gaps lie between secondary structure segments or inside them (Figure 9). Based on the fact that the gaps do not disrupt a longer continuous secondary structure segment, one could argue that the gaps lie between ordered secondary structure segments. But given these specific alignments, there is not enough evidence to support this statement.

In our opinion, there are two distinct quality measures for multiple alignments. The first one is the number of highly conserved columns and the second one is the number of gaps in the alignment. There should be as many conserved columns and as less gaps as possible. These two measures should always be considered together and not isolated since the input parameters can be modified to maximize one of the two parameters. Therefore, a tradeoff between the two criteria has to be found in practice that always depends on the desired results. For example, if one searches for signal peptides, i.e. short and highly conserved sequence stretches, lower gap penalties should be chosen as when comparing the overall similarity of proteins.

Based on our experience, we would choose ClustalW for multiple alignments, because it yielded the best tradeoff between usability and result quality.

Difference between revisions of "Task 2: Alignments"

Latest revision as of 13:38, 1 September 2013

Contents

Sequence Searches

Blast

Psiblast

PDB search

HHblits

Comparison

Multiple sequence alignments

ClustalW

MAFFT

T-Coffee

Expresso

Discussion

Navigation menu

Views

Personal tools

Bioinformatik navigation

MediaWiki navigation

Search

Tools

@@ Line 1: / Line 1: @@
-Sequence alignments are a good start to analyse a new protein sequence. Alignment search tools such as Blast, Psiblast or hhblits can be used to find related protein sequences.
+Sequence alignments are a good start to analyse a new protein sequence. Alignment search tools such as Blast, Psiblast or HHblits can be used to find related protein sequences.
+[[Lab journal task 2]]
 == Sequence Searches ==
-All searches whre done with an increase number of output lines (10000) in the summary hit list and for the reported alignments in order to have all found hits displayed. The analyses were all conducted using all displayed hits. We set no specific evalue cutoff.
-The HFE protein sequence has the Uniprot ID Q30201, NCBI ID 1890180 and PDB ID 1A6Z_A.
+The HFE protein sequence was used to conduct sequence searches with the three different search tools Blast, Psiblast and HHblits.
+CATH and COPS, two protein structure classification databases, were used to evaluate and compare the results of the different tools.
-All Blast, Psiblast and hhblits output files that where analyses where first parse using the perl script [[parse_output.pl]]. For example:
- perl /mnt/home/student/kalemanovm/master_practical/Assignment2_Alignments/scripts/task1/parse_output.pl --out_p /mnt/home/student/betza/task2/blast/res_blast.txt
 CATH is a database that hierarchically classifies all domain structures from PDB. This hierarchy consists of four major levels:
@@ Line 17: / Line 16: @@
 We compared the first level of the CATH hierarchy (Class) to check whether the hits have the same fold class as the query.
 The class of a domain is determined by its secondary structure composition (mostly alpha, mostly beta, mixed alpha/beta or few secondary structures).
-Each domain of a protein is assigned to a CATH class, this means, that one query protein can consist of several fold classes, one for each domain.
+Each domain of a protein is assigned a CATH class, this means, that one query protein can consist of several fold classes, one for each domain.
-The python script [[compareCath.py]] was written to check the overlap of the query protein's fold classes with those of the hits. The two domains of the HFE protein belong to fold class 2 and 3.
+The HFE protein consists of two domains that belong to the CATH fold class 2 and 3.
+Another database for protein structure clustering is COPS. It consists of different structural groups that contain all proteins with a structural similarity above a defined threshold. For example, all proteins in the L30 group have a structural similarity of at least 30% percent. This means that if two proteins are not in the same L30 group, they have less than 30% structural similarity.
-Example call:
+The L30 groups of the hits were determined to check if they are in the same group as the HFE protein or not.
- python /mnt/home/student/betza/scripts/compareCath.py -i /mnt/home/student/betza/task2/blast/res_blast.txt_results -q 1a6zA > /mnt/home/student/betza/task2/blast
- /res_blast.txt/blast_cath
-Another database for protein structure clustering is COPS. It consists of different structural groups that contain all proteins with a structural similarity above athreshold. For example, all proteins in the L30 group have a structural similarity of at least 30% percent. This means that if two proteins are not in the same L30 group, they have less than 30% structural similararity.
-The L30 groups of the hits were determined to chek if they are in the same group as HFE or not. This analyses were also done using the script [[parse_output.pl]]:
- perl /mnt/home/student/kalemanovm/master_practical/Assignment2_Alignments/scripts/task1/parse_output.pl --out_p /mnt/home/student/betza/task2/blast/res_blast.txt
- --query 1a6z_A --sot L30
 ===Blast===
-A Blast search in big_80 was executed using the standard parameter settings:
+The Blast search in the big_80 database yielded 1504 hits.
+<figure id="blast distribution">
- blastall -p blastp -i /mnt/home/student/betza/data/hfe.fasta -d /mnt/project/pracstrucfunc13/data/big/big_80 -o /mnt/home/student/betza/task2/blast/res_blast.txt
+{| align="center"
- -v 10000 -b 10000
+| align="center" | [[File:blast_hist_e.png|center|thumb|450px|'''a)''' E-value distribution of the results of the Blast search against big_80 using standard parameters. The x-axis shows the base 10 logarithm of the e-value and the y-axis the frequency. The better the e-value the smaller the log(e-value).]]
+| align="center" | [[File:blast_hist_i.png|center|thumb|450px|'''b)''' Identity distribution of the results of the Blast search against big_80 using standard parameters. The x-axis shows the relative identity between the query and the hits.]]
+|+ style="caption-side: bottom; text-align: left" |<font size=1.5>'''Figure 1:''' E-value and identity distribution of the Blast search hits.
+|}
+</figure>
+<xr id="blast distribution"/> shows the e-value and identity distributions of all unfiltered Blast search hits.
-Blast finds 1504 hits. There are only a few hits with a very high evalue of 1e-160 but most good hits have an evalue between 1e-40 and 1e-50 (left peak). Most results have a worse evalue above 1e-20. Nearly all hits have an identity of 20% to 40% with the query. The evalue and identity distributions are shown below in the plots.
+There are only a few hits with a very low e-value of 1e-160, but most significant hits have an e-value between 1e-40 and 1e-50 (left small peak in a)). Most results have a worse e-value above 1e-20. The histogram of the sequence identity b) shows that nearly all hits have an identity between 20% to 40% to the query.
+The Blast results only contain 13 PDB IDs. Nevertheless, the number of shared CATH classes between the reference sequence HFE and the hits were computed as described in the [[Lab journal task 2| lab journal]].
-[[File:blast_hist_e.png|left|frame|Evalue distribution of the results of the Blast search against big_80 using standard parameters. The x-axis shows the natural logarithm of the evalue and the y-axis the frequency. The better the evalue the smaller is log(evalue).]]
-[[File:blast_hist_i.png|center|frame|Identity distribution of the results of the Blast search against big_80 using standard parameters. The x-axis shows the relative identity between the query and the hits.]]
+<css>
-The Blast results only contain 13 PDB IDs. Nevertheless, the number of shared CATH classes betwen HFE and the hits were computed as described above.
+table.colBasic2  {
-Of 13 pdb hits, only 3 share all domain fild classes with HFE, 6 have one domain fold class in commom and 4 even none.
+margin-left: auto;
+margin-right: auto;
+border: 2px solid black;
+border-collapse:collapse;
+width: 40%;
+}
+.colBasic2 th,td {
+padding: 3px;
+border: 2px solid black;
+}
+.colBasic2 td {
+text-align:left;
+}
+.colBasic2 tr th   { background-color:#efefef; color: black;}
+.colBasic2 tr:first-child th { background-color:#adceff; color:black;}
+</css>
+<figtable id="blast cath">
-{| border="1"
+{|class="colBasic2" align="right"
-|+ Histogram of the number of shared CATH classes between HFE and the Blast hits.
-! number of same fold classes !! frequency
+! number of same fold classes || frequency
 |-
 ! 0
@@ Line 57: / Line 69: @@
 ! 2
 | 3
+|+ style="caption-side: bottom; text-align: left" |<font size=2>'''Table 1:''' Histogram of the number of shared CATH classes between HFE and the Blast hits.
 |}
+</figtable>
+<xr id="blast cath"/> shows, that only 3 of the 13 PDB hits share all domain fold classes with HFE, 6 of the hits have one domain fold class in commom with the reference and 4 even none. Besides, only 6 of the 13 pdb hits are in the same COPS L30 group as the HFE protein.
-Only 6 of the 13 pdb hits are in the same COPS L30 group as the HFE protein.
+Those numbers indicate, that some of the Blast hits are not really related to HFE and thus, false positives.
-All results indicate, that some of the Blast hits are not really related to HFE and thus false positives.
 ===Psiblast===
-Four different psiblst runsagaisnt big_80 where done using all four combinaitons of:
+Four different Psiblast runs against big_80 were executed using the following four parameter combinations:
-*2 iterations
+*2 iterations and default e-value cutoff 0.002 (j=2, h=2)
-*10 iterations
+*2 iterations and e-value cutoff 10E-10  (j=2, h=10E-10)
-*default E-value cutoff (0.002)
+*10 iterations and default e-value cutoff 0.002 (j=10, h=2)
-*E-value cutoff 10E-10
+*10 iterations and e-value cutoff 10E-10 (j=10, h=10E-10)
+The e-value and sequence identity distributions of the hits from the four runs were compared to identify the best parameter combination.
-Example call for Psiblast with 10 iterations and evalue cutoff of 10E-10:<br>
- blastpgp -i /mnt/home/student/betza/data/hfe.fasta -d /mnt/project/pracstrucfunc13/data/big/big_80 -j 10 -h 10E-10 -v 10000 -b 10000 -o /mnt/home/student/betza
- /task2/psiblast/new/j1h1/j1hi.txt -Q /mnt/home/student/betza/task2/psiblast/new/j1h1/j1h1.pssm -C /mnt/home/student/betza/task2/psiblast/new/j1h1/j1h1.chk
+<figure id="psiblast boxplot">
-The perl scrip devide_psiblast_out.pl was then used to divide the different psiblast iterations in order to be able to analyse the results of the last iteration alone.
+{| align="center"
- perl /mnt/home/student/kalemanovm/master_practical/Assignment2_Alignments/scripts/task1/devide_psiblast_out.pl /mnt/home/student/betza/task2/psiblast/new/j1h1
+|align="center" |  [[File:psi_box_e.png|center|thumb|450px|'''a)''' Boxplot of the log(e-value) of the four different Psiblast searches against big_80: 2 iterations with default e-value cutoff (0.002) and 10E-10 and 10 iterations with default e-value cutoff (0.002) and 10E-10. The thick black line inside the box is the median and the upper and lower end of the box the 25% and 75% quantils. The single points above the whiskers are outliers.]]
- /j1hi.txt
+|align="center" |  [[File:psi_box_i.png|center|thumb|450px|'''b)''' Boxplot of the identity of the results of the four different Psiblast searches against big_80: 2 iterations with default e-value cutoff (0.002) and 10E-10 and 10 iterations with default e-value cutoff (0.002) and 10E-10. The thick black line inside the box is the median and the upper and lower end of the box the 25% and 75% quantils. The single points above and below the whiskers are outliers.]]
+|+ style="caption-side: bottom; text-align: left" |<font size=1.5>'''Figure 2:''' E-value and percent sequence identity distribution of the results from the four Psiblast runs against big_80.
+|}
+</figure>
+<figtable id="psiblast summary">
-The evalue an identity distributions of the four runs are plotted below.
+{|class="colBasic2" align="right"
+!colspan="2" |  || colspan="4" |Psiblast parameters
-[[File:psi_hist_4.png|left|frame|Histogram of the log(evalue) distribution of the results of the four different Psiblast searches against big_80: 2 iterations with default evalue cutoff (0.002) and 10E-10 and 10 iterations with default evalue cutoff(0.002) and 10E-10. The x-axis shows the natural logarithm of the evalue and the y-axis the frequency. The better the evalue the smaller is log(evalue).]]
+|-
+!colspan="2" | ||  j=2, h=0.002 || j=2, h=10E-10 || j=10, h=0.002 || j=10, h=10E-10
-[[File:psi_box_e.png|center|frame|Boxplot of the log(evalue) of the four different Psiblast searches against big_80: 2 iterations with default evalue cutoff (0.002) and 10E-10 and 10 iterations with default evalue cutoff(0.002) and 10E-10. The thick black line inside the box is the median and the upper and lower end of the box the 25% and 75% quantils. The single points above the whiskers are outliers.]]
-[[File:psi_hist_i.png|left|frame|Histogram of the identity distribution of the results of the four different Psiblast search against big_80: 2 iterations with default evalue cutoff (0.002) and 10E-10 and 10 iterations with default evalue cutoff(0.002) and 10E-10. The x-axis shows the relative identity between the query and the hits.]]
-[[File:psi_box_i.png|center|frame|Boxplot of the identity of the results of the four different Psiblast search against big_80: 2 iterations with default evalue cutoff (0.002) and 10E-10 and 10 iterations with default evalue cutoff(0.002) and 10E-10. The thick black line inside the box is the median and the upper and lower end of the box the 25% and 75% quantils. The single points above and below the whiskers are outliers.]]
-The boxplot of the evalue distribution shows, that the default parameter settings (j=2 and h=0.002) can be considered as the best. The median of log(evalue) is higher than in all the other 3 psiblast runs. This indicates that the percentage of hits with a good evalue is higher in the default psiblast run than in the other three.
-The results of the CATH and COPS classification analysis are listed below:
-{| border="1"
-|+ Comparison of the results of the four Psiblast searches against big_80.
-! psiblast parameters !! j=2, h=0.002 !! j=2, h=10E-10 !! j=10, h=0.002 !! j=10, h=10E-10
 |-
-! # hits
+! colspan="2" | # hits
 | 2564 || 3163 || 2601 || 3343
 |-
-! # PDB hits
+! colspan="2" | # PDB hits
-| 66 || 72 || 56 || 66
+|  66 || 72 || 56 || 66
 |-
+! rowspan="3" | CATH classes in common with HFE
-! 0 CATH class
+!  0
 | 21 || 21 || 20 || 21
 |-
-! 1 CATH class
+! 1
 | 42 || 48 || 33 || 42
 |-
+!2
-! 2 CATH classes
 | 3 || 3 || 3 || 3
 |-
-!same L30 group
+!  colspan="2" | # hits within same L30 group
 | 6 || 6 || 6 || 6
+|+ style="caption-side: bottom; text-align: left" |<font size=2>'''Table 2:''' Comparison of the results of the four Psiblast searches against big_80.
-|-
 |}
+</figtable>
+<xr id="psiblast summary"/> summarises the results of the different Psiblast runs. The run with default parameter settings (j=2 and h=0.002) found the least hits.
-There are 1095 proteins in the same COPS L30 group as the query HFE. Since the big_80 datatbase does not contain all proteins from pdb, the quality of the program cannot be assesed using the COPS L30 group of HFE. Nevertheless, of the found pdb hits only a few belong to the same L30 and share the same CATH fold classes with HFE.
+Nevertheless, the boxplot of the e-value distribution (<xr id="psiblast boxplot"/> a)) shows, that the default parameters can be considered the best: the median of the log(e-value) is smaller than the median of the other 3 psiblast runs. This indicates that the percentage of hits with a low and thus significant e-value is higher in the default psiblast run than in the other three.
+Besides, the default parameters also led to a higher percent sequence identity as indicated by the median of the cyan box in <xr id="psiblast boxplot"/> b).
+The number of PDB hits varies among the different runs, but all found 3 proteins with the same two CATH classes as the query.
+There are 1095 proteins in the same COPS L30 group as the query HFE. Since the big_80 database does not contain all proteins from the PDB, the quality of the program cannot be assessed using the COPS L30 group of HFE. Nevertheless, from all found PDB hits, only 6 proteins belonged to the same L30 group.
+====PDB search====
-To get more pdb ids, we also did 4 psiblast searches against /mnt/home/rost/kloppmann/data/blast_db/pdb_seqres using the profiles created with big_80. Those four results were compared as well.
+To get more PDB ids, we also did 4 Psiblast searches against the PDB seqres database using the profiles created with big_80. Those four result lists were compared as well.
-[[File:pdb_psi_hist_e.png|left|frame|Histogram of the log(evalue) distribution of the results of the four Psiblast searches against pdb_seqres using the checkpoint files from the runs before. (2 iterations with default evalue cutoff (0.002) and 10E-10 and 10 iterations with default evalue cutoff(0.002) and 10E-10). The x-axis shows the natural logarithm of the evalue and the y-axis the frequency. The better the evalue the smaller is log(evalue).]]
-[[File:pdb_psi_box_e.png|center|frame|Boxplot of the log(evalue) of the four different Psiblast searches against pdb_seqres using the checkpoint files from the runs before. (2 iterations with default evalue cutoff (0.002) and 10E-10 and 10 iterations with default evalue cutoff(0.002) and 10E-10). The thick black line inside the box is the median and the upper and lower end of the box the 25% and 75% quantils. The single points below the whiskers are outliers.]]
+<figure id="psiblast pdb">
+{| align="center"
+|align="center" |  [[File:pdb_psi_box_e.png|center|thumb|450px|''' b)''' Boxplot of the log(e-value) of the four different Psiblast searches against PDB_seqres using the checkpoint files from the earlier runs. (2 iterations with default e-value cutoff (0.002) and 10E-10 and 10 iterations with default e-value cutoff(0.002) and 10E-10). The thick black line inside the box is the median and the upper and lower end of the box the 25% and 75% quantils. The single points below the whiskers are outliers.]]
+|align="center" | [[File:pdb_psi_box_i.png|center|thumb|450px|Boxplot of the identity of the results of the four different Psiblast searches against pdb_seqres using the checkpoint files from the runs before. (2 iterations with default e-value cutoff (0.002) and 10E-10 and 10 iterations with default e-value cutoff(0.002) and 10E-10). The thick black line inside the box is the median and the upper and lower end of the box the 25% and 75% quantils. The single points above and below the whiskers are outliers.]]
+|+ style="caption-side: bottom; text-align: left" |<font size=1.5>'''Figure 3:''' E-value and sequence identity distributions of the hits from the four Psiblast runs against PDB_seqres.
+|}
+</figure>
+<figtable id="psiblast summary">
-[[File:pdb_psi_hist_i.png|left|frame|Histogram of the identity distribution of the results of the four different Psiblast searches against pdb_seqres using the checkpoint files from the runs before. (2 iterations with default evalue cutoff (0.002) and 10E-10 and 10 iterations with default evalue cutoff(0.002) and 10E-10). The x-axis shows the relative identity between the query and the hits.]]
+{|class="colBasic2" align="center"
+!colspan="2" |  || colspan="4" |Psiblast parameters
-[[File:pdb_psi_box_i.png|center|frame|Boxplot of the identity of the results of the four different Psiblast searches against pdb_seqres using the checkpoint files from the runs before. (2 iterations with default evalue cutoff (0.002) and 10E-10 and 10 iterations with default evalue cutoff(0.002) and 10E-10). The thick black line inside the box is the median and the upper and lower end of the box the 25% and 75% quantils. The single points above and below the whiskers are outliers.]]
-{| border="1"
-|+ Comparison of the results of the four Psiblast searches against pdb_seqres.
-! psiblast checkfiles !! j=2, h=0.002 !! j=2, h=10E-10 !! j=10, h=0.002 !! j=10, h=10E-10
 |-
+! colspan="2" | ||  j=2, h=0.002 || j=2, h=10E-10 || j=10, h=0.002 || j=10, h=10E-10
-! # hits
+|-
+!colspan="2" |  # hits
 | 4307 || 4814 || 3874 || 4058
 |-
+! rowspan="3" | CATH classes in common with HFE
-! 0 CATH class
+! 0
 | 561 || 586 || 529 || 539
 |-
-! 1 CATH class
+! 1
 | 2798 || 3280 || 2397 || 2571
 |-
-! 2 CATH classes
+! 2
 | 948 || 948 || 948 || 948
 |-
-!same L30 group
+!colspan="2" |  # hits within same L30 group
 | 1047 || 1053 || 1045 || 1052
 |-
+|+ style="caption-side: bottom; text-align: left" |<font size=2>'''Table 3:''' Comparison of the results of the four Psiblast searches against PDB_seqres.
 |}
+</figtable>
+The e-value and sequence identity distributions in <xr id="psiblast pdb"/> do not show a significant difference in the results, because the medians of all four e-value distributions are nearly the same. The last run with the parameter combination of j=10 and h=10E-10 has a higher percentage of low e-value hits, because the 25% quantil is below those of the other distributions. But the run with j=2 and h=10E-10 lead to more hits with a higher sequence identity, because the median slightly lies above the other medians. Apart from the query sequence, there were no PDB hits with a sequence identity above 40%.
-proteins are in the same COPS L30 group as HFE. The differnent psiblast searches found nearly all proteins that are in the same L30 group as HFE. Therefore,the
-decision which of the parameter settings leads to the best results cannot be based ont he pdb hits, since the differences are minimal.Nevertheless, j=10, h=0.002 had the least hits, but the highest percentage of hits in the same L30 group as the query.
+Evaluating the CATH and SCOP annotations of the found hits revealed that all four searches found 948 hits with the same 2 CATH classes as reference. Besides, all runs could find nearly all 1095 proteins that are in the same COPS L30 group as HFE. Therefore, the decision which of the parameter settings leads to the best results cannot be based on the PDB hits, since the differences are minimal. Nevertheless, j=10, h=0.002 had the least hits, but the highest percentage of hits in the same L30 group as the query. It is thus best suited to find a high number of significant hits.
+<figure id="psiblast venn">
-We also checked how many ot the hits overlapped in the four different runs. The results are shown in the Venn diagramm below:
+{| align="center"
+|align="center" |  [[File:Psiblast_venn_comp.png|center|thumb|400px]]
+|+ style="caption-side: bottom; text-align: left" |<font size=1.5>'''Figure 4:''' Venn diagram of the overlap between all hits from the four different searches. J1H1: 10 iterations and e-value cutoff 10E-10, J1H2: 10 iterations and e-value cutoff 2E-3, J2H1: 2 iterations and e-value cutoff 10E-10, J2H2: 2 iterations and e-value cutoff 2E-3.
+|}
+</figure>
+We also computed the hits overlap of the four different Psiblast searches in big_80, to see how comparable the found hits are. The overlap is illustrated in the Venn diagram in <xr id="psiblast venn"/>.
-TODO Venn diagramm !!!!
-Additionally, the 100 best hits of each of the runs were compared to the 100 best hits of all the other 3 runs. They were the same with all four parameter combinations. Thus, the parameter settings only matter if a high number of hits is needed, but not of onyl the best hits are used in further analyses.
+Additionally, we computed the overlap for the 100 best hits (based on e-value) of each of the runs (data not shown). They were exactly the same in all four runs. Thus, the parameter settings only matter if a high number of hits is needed, but not if onyl the best hits of a search are used in further analyses.
 ===HHblits===
-hhblits -i /mnt/home/student/betza/data/hfe.fasta -d /mnt/project/rost_db/data/hhblits/uniprot20_02Sep11 -o /mnt/home/student/betza/task2/hhblits/hfe.hhr -oa3m /mnt/home/student/betza/task2/hhblits/hfe.a3m -oalis /mnt/home/student/betza/task2/hhblits/hfe -ohhm /mnt/home/student/betza/task2/hfe.hhm -Z 10000 -B 10000
+We executed HHblits with default parameters and searched the uniprot20_02Sep11 database. HHblits is an alignment search tool that searches in a database of clustered sequences. A profile is generated from the query that is compared with the profile of each cluster. For each cluster, not only the representative, but also all other cluster members are reported in the results list.
-The output files were analysed using the script parse_output.pl.
+Thus, HHblits finds much more sequences than Blast and Psiblast. We found 2861 clusters and 37908 sequences in total. We used all sequences for the analysis, because one of HHblits' most important features is the high number of good results.
-example:
-perl /mnt/home/student/kalemanovm/master_practical/Assignment2_Alignments/scripts/task1/parse_output.pl --out_p /mnt/home/student/betza/task2/blast/res_blast.txt --query 1a6z_A --sot L30 --out_h /mnt/home/student/betza/task2/hhblits/
+<figure id="hhblits distribution">
-[[File:hhblits_hist_e.png|left|frame|Evalue distribution of the results of the HHblits search with default parameters using the using the uniprot20_02Sep11 database. The x-axis shows the natural logarithm of the evalue and the y-axis the frequency. The better the evalue the smaller is log(evalue).]]
+{| align="center"
+| align="center" | [[File:hhblits_hist_e.png|center|thumb|450px|'''a)''' E-value distribution of the results of the HHblits search with default parameters using the uniprot20_02Sep11 database. The x-axis shows the base 10 logarithm of the e-value and the y-axis the frequency. A smaller log10(e-value) corresponds to a better results.]]
+| align="center" | [[File:hhblits_hist_i.png|center|thumb|450px|'''b)''' Identity distribution of the results of the HHblits search in the uniprot20_02Sep11 database using default parameters. The x-axis shows the relative identity between the query and the hits.]]
+|+ style="caption-side: bottom; text-align: left" |<font size=1.5>'''Figure 5:''' E-value and identity distributions of the HHblits search hits.
+|}
+</figure>
+<xr id ="hhblits distribution"/> a) shows the distribution of the log10 of the e-values of the HHblits hits. It is remarkable that HHblits finds so many hits with a very low e-value of 1E-130. One reason is, that the best hit belongs to a very large cluster. Nevertheless, in comparison to the Blast e-value distribution ( <xr id="blast distribution"/> a)), the percentage of hits with a high e-value is much lower in the HHblits results.
-[[File:hhblits_hist_i.png|center|frame|Identity distribution of the results of the HHblits search in the uniprot20_02Sep11 database using default parameters. The x-axis shows the relative identity between the query and the hits.]]
+There are several hits with a very low sequence identity to the query (see <xr id="hhblits distribution"/> b)), but the distribution contains two peaks around 20% and 40% identity that correspond to distant and more closely related hits.
-It is remarkable that HHblits finds so many hits with a very low evalue. This can maybe explained by the fact, that hhblits found several clusters with a good evalue. Since used all proteins in the clusters for the analyses and not only the cluster representative the number of hits is clearly higher than the number of hits from the Blast and Psiblast searches.
-===Comparison===
+==Comparison==
+<figure id="comparison">
-[[File:compDefault_e.png|left|frame|TODO evalue]]
+{| align="center"
+| align="center" | [[File:compDefault_e.png|center|thumb|450px|'''a)''' E-value distribution of the Blast, Psiblast and HHblits hits with default parameter searches.]]
+| align="center" | [[File:compDefault_i.png|center|thumb|450px|'''b)''' Identity distribution of the Blast, Psiblast and HHblits hits default parameters searches]]
+|+ style="caption-side: bottom; text-align: left" |<font size=1.5>'''Figure 6:''' E-value and identity distribution of the Blast, Psiblast and HHblits search hits.
+|}
+</figure>
+We used the results of the Psiblast run with 2 iterations and an e-value cutoff of 2e-3, because we also used default parameters for Blast and HHblits and because we also received the best results with this parameters.
-[[File:compDefault_i.png|center|frame|TODO idenity]]
+We compared the e-value and identity distributions. <xr id ="comparison"/> a) clearly shows that Blast has the lowest e-values, but the percentage of low e-value hits is higher in the HHblits results. The median sequence identity has the highest value among the Blast hits.
+<figure id="comparison venn">
+{| align="center"
+|align="center" |  [[File:Method_comp_venn.png|center|thumb|400px]]
+|+ style="caption-side: bottom; text-align: left" |<font size=1.5>'''Figure 7:''' Venn diagram of the overlap between all hits from the Blast, Psiblast and HHblits runs with default parameters.
+|}
+</figure>
+The overlap between all hits and the top 100 hits are shown in <xr id ="comparison venn"/>.
+<figtable id="best hits">
-TODO: Venn diagramms of hits overlap.
+{|class="colBasic2" align="right"
+! search tool || colspan=6 |rank of hit
+|-
+! Psiblast
+| 1   ||  2     ||  3 ||      4 ||      5    ||   6
+|-
+! Blast
+|  1   ||   3 ||     2   ||  52 ||   109  ||   81
+|-
+! HHblits
+| 2  ||   not found   ||    not found  ||     not found  ||     not found      || not found
+|+ style="caption-side: bottom; text-align: left" |<font size=2>'''Table 4:'''  Rank of the top 6 Psiblast hits in the Blast and HHblits results.
+|}
+</figtable>
+We tried to identify the 6 best Psiblast hits with the lowest e-values in the results of the Blast and HHblits searches, see <xr id ="best hits"/>-
+The top 6 Psiblast hits are all present in the Blast results and the first 3 hits are the same. This confirms that Blast and Psiblast are very similar.
+HHblits, on the contrary, has definitely not the same results as Blast and Psiblast. From the top 6 Psiblast hits, only the first one is present in the second best cluster. This proves, that HHblits finds more distant related hits, which is also indicated by the low median sequence identity of the HHblits hits.
+The best 5 hits from all three search tools belong to the MHC class I family and are thus correct hits.
 == Multiple sequence alignments==
+{| align="center" style="border-collapse:collapse"
-The following data sets were chosen for the Multiple Sequence Alignments:
+| style="background-color:#adceff; text-align:center; font-weight:bold" colspan="2" | < 30% sequence identity || colspan="2" style="background-color:#adceff; text-align:center; font-weight:bold" | > 60% sequence identity || colspan="4" style="background-color:#adceff; text-align:center; font-weight:bold" | whole range of sequence identity
+|-
+| style="background-color:#adceff; text-align:center; font-weight:bold" | AC Number || style="background-color:#adceff; text-align:center; font-weight:bold" | Sequence Identity || style="background-color:#adceff; text-align:center; font-weight:bold" | AC Number || style="background-color:#adceff; text-align:center; font-weight:bold" | Sequence Identity || style="background-color:#adceff; text-align:center; font-weight:bold" | AC Number || style="background-color:#adceff; text-align:center; font-weight:bold" | Sequence Identity || style="background-color:#adceff; text-align:center; font-weight:bold" | AC Number || style="background-color:#adceff; text-align:center; font-weight:bold" | Sequence Identity
+|-
+| F6JYA9 || 0.29 || B4DDZ1 || 0.83 || B4DDZ1 || 0.83 || F6WCX4 || 0.34
+|-               ||                ||                ||
+| D9J389 || 0.28 || Q5EEZ1 || 0.76 || Q5EEZ1 || 0.76 || F6JYA9 || 0.29
+|-               ||                ||                ||
+| D5MSB3 || 0.27 || G3THV5 || 0.75 || G1MBW1 || 0.73 || D9J389 || 0.28
+|-               ||                ||                ||
+| B3FRK2 || 0.20 || G1MBW1 || 0.73 || H0VAR7 || 0.72 || D5MSB3 || 0.27
+|-               ||                ||                ||
+| 3ov6_A || 0.22 || H0VAR7 || 0.72 || F1PX48 || 0.71 || 2zok_E || 0.24
+|-               ||                ||                ||
+| 1p7k_L || 0.21 || F1PX48 || 0.71 || G1T7D7 || 0.70 || H0Y1D0 || 0.20
+|-               ||                ||                ||
+| H0Y1D0 || 0.20 || G1T7D7 || 0.70 || G5BQE5 || 0.67 || Q8SNJ4 || 0.22
+|-               ||                ||                ||
+| B3FRK3 || 0.18 || O35799 || 0.68 || Q95IT9 || 0.38 || B3FRK3 || 0.18
+|-               ||                ||                ||
+| Q8HWL2 || 0.17 || G5BQE5 ||0.67  || 2qrt_A || 0.38 || Q8HWL2 || 0.17
+|-                                 ||
+|  ||  ||  ||  || Q8HX83 || 0.36    || ||
+|+ style="caption-side: bottom; text-align: left" |<font size=2>'''Table 5:''' List of all sequences that were used in the multiple sequence alignments.
+|}
+In order to assess the difference of multiple sequence alignments between close and related homologs, three different groups of sequences were selected. One containing only sequences with a sequence identity above 60% to human hfe, one with sequences below 30% identity and one with sequences covering the whole range of sequence identity. The selected sequences are listed in Table 5. <br />
-Below 30% SeqID:
+In order to ensure that also sequences with known structures are included in the alignments, the sequences from the following PDB structures were selected:
+*1p7k_L
+*3ov6_A
+*2qrt_A
+*2zok_E
-{| class="wikitable"
+===ClustalW===
+<figtable id="clustalW">
+{| style="border-style: none ; border-collapse:collapse; border-width: 0px"
+| [[File:Below30_cw.png|center|thumb|500px| Below 30% identity.]] || [[File:Above60_cw.png|center|thumb|500px| Above 60% identity.]] || [[File:wholerange_cw.png|center|thumb|500px| Whole range of identity.]]
 |-
+|+ style="caption-side: bottom; text-align: left" |<font size=2>'''Figure 8:''' ClustalW alignments for all different groups of sequence identity are shown.
-| AC Number || Sequence Identity
 |-
+|}
-| F6JYA9 || 0.29
+</figtable>
+<figtable id="clustalWSS">
+{| style="border-style: none ; border-collapse:collapse; border-width: 0px"
+| [[File:Clustalw_below30_ss.png|center|thumb|500px| ClustalW alignment of sequences from the below 30% identity group.]] || [[File:Clustalw_above60_ss.png|center|thumb|500px| ClustalW alignment of sequences from the above 60% identity group.]] || [[File:Clustalw_wholerange_ss.png|center|thumb|500px| ClustalW alignment of sequences from the whole range identity group.]]
 |-
+|+ style="caption-side: bottom; text-align: left" |<font size=2>'''Figure 9:''' The multiple alignments are displayed with the secondary structure annotation of Q30201 from Uniprot shown in red (helix), blue (sheet) and green (loop).
-| D9J389 || 0.28
-|-
+|}
+</figtable>
-| D5MSB3 || 0.27
-|-
+===MAFFT===
-| B3FRK2 || 0.20
-|-
-| 3ov6_A || 0.22
+<figtable id="mafft">
-|-
+{|
-| 1p7k_L || 0.21
+| [[File:Below30_mafft.png|center|thumb|500px| Below 30% identity.]] || [[File:Above60_cw.png|center|thumb|500px| Above 60% identity.]] || [[File:Whole_mafft.png|center|thumb|500px| Whole range of identity.]]
-|-
-| H0Y1D0 || 0.20
 |-
+|+ style="caption-side: bottom; text-align: left" |<font size=2>'''Figure 10:''' MAFFT alignments for all different groups of sequence identity are shown.
-| B3FRK3 || 0.18
 |-
-| Q8HWL2 || 0.17
 |}
+</figtable>
-Above 60% SeqID:
+===T-Coffee===
-{| class="wikitable"
+<figtable id="tcoffee">
+{|
+| [[File:Below30_tcof.png|center|thumb|500px| Below 30% identity.]] || [[File:Above60_tcof.png|center|thumb|500px| Above 60% identity.]] || [[File:Whole_tcof.png|center|thumb|500px| Whole range of identity.]]
 |-
+|+ style="caption-side: bottom; text-align: left" |<font size=2>'''Figure 11:''' T-Coffee alignments for all different groups of sequence identity are shown.
-| AC Number || Sequence Identity
 |-
-| B4DDZ1 || 0.83
-|-
-| Q5EEZ1 || 0.76
-|-
-| G3THV5 || 0.75
-|-
-| G1MBW1 || 0.73
-|-
-| H0VAR7 || 0.72
-|-
-| F1PX48 || 0.71
-|-
-| G1T7D7 || 0.70
-|-
-| O35799 || 0.68
-|-
-G5BQE5	0.67
 |}
+</figtable>
+===Expresso===
-Whole range of SeqID:
+<figtable id="tcoffee">
-{| class="wikitable"
+{|
+| [[File:Below30_ex.png|center|thumb|500px| Below 30% identity.]] || [[File:Above60_ex.png|center|thumb|500px| Above 60% identity.]] || [[File:Wholerange_ex.png|center|thumb|500px| Whole range of identity.]]
 |-
+|+ style="caption-side: bottom; text-align: left" |<font size=2>'''Figure 12:''' Expresso alignments for all different groups of sequence identity are shown.
-| AC Number || Sequence Identity
 |-
-| B4DDZ1 || 0.83
-|-
-| Q5EEZ1 || 0.76
-|-
-| G1MBW1 || 0.73
-|-
-| H0VAR7 || 0.72
-|-
-| F1PX48 || 0.71
-|-
-| G1T7D7 || 0.70
-|-
-| G5BQE5 || 0.67
-|-
-| Q95IT9 || 0.38
-|-
-| 2qrt-A || 0.38
-|-
-| Q8HX83 || 0.36
-|-
-| F6WCX4 || 0.34
-|-
-| F6JYA9 || 0.29
-|-
-| D9J389 || 0.28
-|-
-| D5MSB3 || 0.27
-|-
-| 2zok_E || 0.24
-|-
-| H0Y1D0 || 0.20
-|-
-| Q8SNJ4 || 0.22
-|-
-| B3FRK3 || 0.18
-|-
-| Q8HWL2 || 0.17
 |}
+</figtable>
-pdb:
-p7k_L 0.21<br>
-ov6_A 0.22 <br>
-qrt_A 0.38<br>
-zok_E 0.24<br>
+===Discussion===
-===ClustalW===
-Below 30% SeqID: [[ClustalW_Below_30]]<br>
-Above 60% SeqID: [[ClustalW_Above_60]]<br>
-Whole range of SeqID: [[ClustalW_whole]]
+In the above 60% sequence identity group, the residues and the gaps are well conserved, especially in the first third of the sequence. In the other two thirds, the sequence conservation drops slightly, but is still at a generally high level.
-===MAFFT===
+In the below 30% identity group, the number of gaps inside the sequences is not much higher than in the above 60% group, but the overall residue conservation is significantly lower. Nevertheless, there are some very well conserved residues in this group that might be functionally important. In the mixed sequence identity group, the gaps are not as well conserved as in the other two groups and the conservation is even lower than in the below 30% identity group. But this effect is probably strengthened by the higher amount of sequences present. Another peculiarity is that the inclusion of 3D information with Expresso (figure 12), yields an alignment with a lot more gaps, but more conserved columns than the other programs. Although this could be a desirable property in an alignment, this behavior can be reproduced with the other programs by modifying the gap penalty parameters.
+The different alignment programs yield comparable results, i.e. none of them is considerably better or worse than the other two. Nevertheless, some differences can be observed. The first notable one is that MAFFT yields a lot less conserved columns for the low sequence identity group than the other programs. Also, the positioning of consecutive gap columns varies between the programs.
-Below 30%:[[MAFFT_30]]<br>
-Above 60%: [[MAFFT_60]]<br>
-Whole range: [[MAFFT_whole]]
+Since there are only one to two bigger gaps in the alignments and the secondary structure segments are often not very long, it is hard to tell if the gaps lie between secondary structure segments or inside them (Figure 9). Based on the fact that the gaps do not disrupt a longer continuous secondary structure segment, one could argue that the gaps lie between ordered secondary structure segments. But given these specific alignments, there is not enough evidence to support this statement.
-===T-Coffee===
+In our opinion, there are two distinct quality measures for multiple alignments. The first one is the number of highly conserved columns and the second one is the number of gaps in the alignment. There should be as many conserved columns and as less gaps as possible. These two measures should always be considered together and not isolated since the input parameters can be modified to maximize one of the two parameters. Therefore, a tradeoff between the two criteria has to be found in practice that always depends on the desired results. For example, if one searches for signal peptides, i.e. short and highly conserved sequence stretches, lower gap penalties should be chosen as when comparing the overall similarity of proteins.
-Below 30%: [[Tcoffee_30]]<br>
-Above 60%:[[Tcoffee_60]]<br>
-Whole range: [[Tcoffee_whole]]
+Based on our experience, we would choose ClustalW for multiple alignments, because it yielded the best tradeoff between usability and result quality.
-===Comparison===