![[Feedback]](/icons/banner/feedbackACT.gif){kind=icon}
![[For Subscribers]](/icons/banner/subscriberACT.gif){kind=icon}
![[Archive]](/icons/banner/archiveACT.gif){kind=icon}
![[Search]](/icons/banner/searchACT.gif){kind=icon}
![[Contents]](/icons/banner/tocACT.gif){kind=icon}
Previous Article | Next Article 
Journal of Clinical Microbiology, January 2002, p. 182-192, Vol. 40, No. 1
0095-1137/02/$04.00+0 DOI: 10.1128/JCM.40.1.182-192.2002
Copyright © 2002, American Society for Microbiology. All Rights Reserved.
Neurovirosis Division, Virology Department, National Institute for Infectious Disease, "Dr. Carlos G. Malbrán," Buenos Aires, Argentina,1 Service of Diagnostic Microbiology, Centro Nacional de Microbiología, Instituto de Salud Carlos III, Majadahonda, Madrid, Spain2
Received 20 February 2001/ Returned for modification 19 August 2001/ Accepted 28 October 2001
 |
ABSTRACT |
|---|
 Top Abstract IntroductionMaterials and MethodsResultsDiscussionReferences |
|---|
|
INTRODUCTION |
|---|
|
TopAbstractIntroductionMaterials and MethodsResultsDiscussionReferences |
|---|
Because of these problems, several methods were developed for the molecular characterization of the genus (1, 5–7, 17). The coupling of reverse transcription and amplification of the enteroviral RNA by PCR, followed by direct sequencing of the amplified products, is the general approach. Some of these methods analyzed the 5' noncoding, the VP2, or the 3D region of the EVs genome; but the sequence did not always correlate with the corresponding serotype (1, 10, 18). The VP1 region was the most suitable target due to the high correlation between serotypes and sequences and the availability of a large database of EVs sequences (17). Two methods were developed for the typing of EV strains by partial sequencing of VP1 from cell cultures (17) or clinical samples (5). These two approaches were based on the sequencing of amplified products and comparison of the sequences with the VP1 sequences in a database of EVs reference strains by pairwise local alignment (17). However, several problems were found for only a handful of field strains, suggesting that this method of analysis could be inconsistent (5).
Although methods for the sequencing of VP1 may vary, there must be widespread agreement as to the method of analysis so that it can be universally applied. Several methods of analysis are quick and standardized: pairwise sequence alignment, multiple-sequence alignment, and phylogenetic reconstruction. Although there are several choices, users generally apply the available methods instead of first performing an evaluation of them. Nevertheless, the election of the analysis method should be an essential part of the design process, since each method of analysis is based on different assumptions.
Previously proposed sequence analysis methods were based on the similarity index of the pairwise local alignment as a measure of identity of the strain sequence (17). The Gap program (Genetics Computer Group [GCG] software package) uses the algorithm of Needleman and Wunsch (16). The same algorithm is implemented by other commercial packages (e.g., the MegAlign program of the DNASTAR package [3]).
Multiple-sequence alignment is an extension of pairwise sequence alignment but is considered a more thorough method of analysis than pairwise sequence alignment. Two programs, a freeware version of the Clustal W program (version 1.81) (23) and the Pileup program, another part of the GCG package, were evaluated.
Phylogenetic reconstruction of gene trees from sequence data is considered the "gold standard" method of molecular analysis. The best results are obtained when the proper method is chosen and the quality control measures are followed. However, there are several available methods based on distinct assumptions.
In the present study, we compared the efficiencies of analysis methods over a range of values for the various parameters. The methods are simple enough to be performed in virology reference laboratories. A completely characterized database of EV strain sequences from different origins was used for comparison.
|
MATERIALS AND METHODS |
|---|
|
TopAbstractIntroductionMaterials and MethodsResultsDiscussionReferences |
|---|
|
Extraction, amplification, and sequencing. Nucleic acids from clinical samples and isolates were extracted as described previously (4). The dried pellet was resuspended in 15 µl of RNase-free sterile water (Sigma Chemical Co., St. Louis, Mo.) and used immediately. A reverse transcription-nested PCR was used to amplify the VP1 genome region, and a 656-bp product (5) which was located between nucleotide 2874 and nucleotide 3529 of echovirus type 9 strain Barty (GenBank accession no. X92886) was obtained.
Cycle sequencing reactions were performed with the Big Dye terminator kit (Perkin-Elmer Applied Biosystems) by using the same primers used for the nested PCR. Raw sequence data were first analyzed with CHROMAS software (version 1.3; C. McCarthy, 1996, Griffith University, Brisbane, Queensland, Australia), and the forward and reverse sequence data for each sample were aligned by using the MegAlign program (DNASTAR Inc. Software, Madison, Wis.) to obtain the final consensus sequence.
EV sequences used. The sequences of 40 original EV strains were used. Another 64 EV sequences were obtained from the GenBank database, and their accession numbers are included in Table 1. VP1 reference strains from a previously described database (19) were included in every analysis method. The 2A region and the primer sequences were excised from the sequences before analysis.
Pairwise alignment analysis. (i) Martinez-Needleman and Wunsch algorithm. The Martinez (11)-Needleman and Wunsch (16) algorithm (M-NW similarity test) is implemented by the MegAlign program of the DNASTAR commercial package. A similarity score between each pair of sequences was obtained manually after sequential pairwise alignment (M-NW similarity test) was performed. The quality of the alignment could not be measured, and the gap opening penalty (GOP) and the gap extension penalty (GEP) values were used as default values.
(ii) Needleman and Wunsch algorithm. The Needleman and Wunsch (16) algorithm (NW similarity test) is implemented by the Gap program of the GCG package (NW-GCG similarity test), and it was performed with both nucleotide and amino acid sequences.
Given that strains of the same serotype rarely have gaps in their alignments, increasing the GOP and GEP parameters will not modify the similarity score for very related sequences but will decrease it for more diverging sequences. Therefore, to allow the separation of the peak of the similarity score for strains of the same serotype from the peak of the similarity score for strains of related serotypes, GOP options were set as default, double, and quadruple (50, 100, and 200, respectively) with GEP options of 3, 10, and 20, respectively. The randomization option was set on 100 to process the quality control of the alignment.
The quality of the alignment was recorded as the difference between the quality of the best alignment (V) and the mean of the Monte Carlo randomized alignments (m) divided by the standard deviation (SD) of the Monte Carlo simulation ([V - m]/SD), as established previously (2). Significance scores above 15 SDs indicated a nearly ideal alignment, while scores above 5 SDs suggested a good alignment and scores below 5 SDs were regarded with caution.
The conditions used to identify the sequences in both nucleotide schemes were the same as those described previously (19). Even though these conditions had been determined for the GCG program, the same criteria were used for the DNASTAR program since both programs are based on the same algorithm.
A sequence identity of 
75% for any EV prototype strain indicated a homologous serotype, provided that the second-highest score was <70% (next closest serotype). A highest score between 70 and 75% or a second-highest score greater than 70% indicated a tentative identification, while a highest score lower than 70% indicated a nonmatching sequence.
Amino acid analysis conditions were set on the basis of a highest score of 90% and a second-highest score of 85%. These cutoff scores were inferred from a frequency distribution graphic of the pairwise sequence analysis scores when the intraserotype peak differentiated (25). The same kind of analysis had been used before for potyvirus classification (25).
Multiple-sequence alignment analysis. Two widely used multiple-sequence alignment programs, the Clustal W program (version 1.81) and the Pileup program, were compared (23). The Pileup program (which is part of the GCG package) is a progressive pairwise alignment program which implements a simplified version of the Feng and Doolittle algorithm (8a). Sequence similarity is used to cluster sequences, and then a dendrogram is constructed by the unweighted pair group method with arithmetic means, which orders the subsequent pairwise alignments based on the method of Needleman and Wunsch (16).
The Clustal W program (the slow, accurate option) is also based on the Feng and Doolitle algorithm but is improved through dynamic assignment of penalties. The Clustal W program varies gap penalties in relation to sites along the sequence and the relative evolutionary distance among sequences (and subsets of sequences), and the dendrogram is constructed by the neighbor joining method. GOP and GEP costs can be specified separately for the pairwise sequence and multiple-sequence alignments, while the delay divergent sequences option delays the inclusion of sequences with less than the specified sequence identity until the most similar sequences are aligned.
The key parameters for both programs are the cost assigned to the opening of the gap (gap cost), the cost for extension of a gap (gap-length cost) (total gap cost = gap cost + [gap-length cost x gap length]), and the delay divergent sequence option, which is available only in the Clustal W program. These parameters were modified to evaluate the performances of the programs. The alignments of the deduced amino acid sequences of all strains were first compared. The methods were then tested by varying the GOP value between 5 and 100 and by using different GEP values (1, 5, and 10 for the Pileup program and 6.66, 20, 100, and 200 for the Clustal W program). A similar analysis was performed for the nucleotide sequences with the Pileup program and a GEP value of 0.3. The Clustal W program delay divergent option was also evaluated. This value was set at 80 since the limit value of the identity score for two strains of the same serotype was estimated to be 75%.
The quality of the alignments should also be assessed in order to evaluate and compare alignment programs. Two different scores were calculated to estimate the quality of the alignment. The sum-of-pairs score (SPS) increases with the number of sequences correctly aligned and is used to determine the abilities of the programs to align some or all sequences within an alignment. The column score (CS) is a binary score which tests the abilities of the programs to align all the sequences correctly. Both scores were previously developed and tested to evaluate multiple-sequence alignment programs (24), but no reference alignment is considered to normalize the corresponding results in the present study.
Finally, the most consistent alignment was chosen to calculate the homology index to identify the serotype. For multiple-sequence alignments, this score for each pair of sequences was calculated as follows:
 |
A homology index higher than 75% assigned a strain to a serotype.
Phylogenetic approach. The most consistent alignment was used as the input for the phylogenetic methods. Phylogenetic relationships were inferred with the DNAPARS program (as a character state method) and the DNADIST or PROTDIST program (as a distance method) and then with the NEIGHBOUR and KITSCH programs (PHYLIP, version 3.57) (8). The PROTPARS, DNAML, and FITCH programs were not evaluated due to the excessive analysis time needed for large numbers of sequences with a common personal computer system. The statistical significance of phylogenies was estimated by using the SEQBOOT program by bootstrap analysis with 100 pseudoreplicate data sets. Only groupings with bootstrap values higher than 65 were considered significant. The consensus tree obtained with CONSENSE was plotted by TREEVIEW (version 1.2; R. Page, 1995, Glasgow, United Kingdom).
Nucleotide sequence accession numbers. The original sequences detected in the present study were submitted to the GenBank sequence database under accession numbers AF252165 to AF252190 and AF290896 to AF290909.
|
RESULTS |
|---|
|
TopAbstractIntroductionMaterials and MethodsResultsDiscussionReferences |
|---|
Pairwise sequence alignment methods. (i) Quality analysis. The quality analysis showed that each program obtained different results. Since the DNASTAR package did not provide an option to perform the quality analysis, the results were compared with these obtained with the GCG package on a pairwise basis to test its quality. One of the most divergent results was observed when isolate R100, identified by seroneutralization as EV type 18 (EV18), was aligned with the reference strain, EV18 Metcalf. The overall length of the Gap alignment was 419 bp, while the MegAlign alignment had an overall length of 531 bp. No gaps were inserted by the Gap program, and the similarity score was 75.3%. On the other hand, the MegAlign program inserted 224 gaps, resulting in a similarity score of 46.3%. The same discordant insertion of gaps was found frequently in other pairs of alignments (data not shown).
The randomization option of the Gap program was used to evaluate the quality of the result by performing a Monte Carlo simulation. The more similar or related the sequences were, the higher the quality of the alignment was (data not shown).
The parameters selected for the analysis of nucleotide sequences were GOP equal to 50 and GEP equal to 3, considering the fact that the lowest-quality scores were 4.99, 2.74, and 1.09 with GOP equal to 50, 100, and 200, respectively (Table 2). By using the selected options, only one alignment resulted in a score lower than 5 (4.99 [1 of 8,216 comparisons; 0.0001%]), while 170 (2.1%) and 946 (11.5%) of the pairwise alignments had scores lower than 5 by using GOP values of 100 and 200, respectively.
|
(ii) Similarity test as method for EV identification. The values obtained by the NW-GCG similarity test of nucleotide sequences from clinical samples are shown in Table 1. Only the EV13 strain and four other field strains (untypeable strains OK85-6388, VA86-6765, CT87-7122, and CT87-7123) resulted in similarities scores <75%, as reported previously (19).
However, the second-highest scores were higher than 70% in more than 60 cases under any condition, resulting in more than 45 different untypeable EV serotypes (Table 2). This was also true for isolate PA88-8412, whose second-highest score was also higher than 75% (the score for strain EV8 Bryson was 75.6%, which was independent of the parameter chosen).
On the other hand, deduced amino acid sequence analysis grouped all tested strains except the untypeable strains with their respective reference strains (including the EV13 field strain). Again, 25 different strains had in second-highest scores higher than 85% (data not shown).
The similarity values obtained by the M-NW similarity test were lower than the values by the NW-GCG similarity test (mean score by the NW-GCG similarity test of 81.5% versus mean score by the M-NW similarity test of 78.8%). A total of 19.2% of the strains, including the untypeable strains, could not be typed because of low similarity scores or because the score attributed the strain to a serotype different from the serotype established by the other methods (Table 3). In conclusion, the M-NW similarity test appears to be less suitable for EV identification.
|
Multiple-sequence alignment methods. (i) Quality analysis. (a) Deduced amino acid sequences. The Pileup program generated multiple-sequence alignments with very different lengths and qualities. The Clustal W program generated alignments with more conserved lengths and qualities for different parameters instead. The performances measured as the overall similarity (SPS/length) and the CSs for the different methods of analysis with different GEP and GOP parameters are shown in Fig. 1a and b.
|
(b) Nucleotide sequences. For the Clustal W program, the alignment of the divergent sequences was delayed with the delay divergent option (which was varied between 30 and 95) until all the remaining sequences were aligned. The effect of this optional parameter in the alignment quality is presented in Fig. 1e, showing that the best alignments were obtained with values higher than 85. Thus, the conditions mentioned above were tested by varying only the values of GOP and GEP. Again, the Clustal W program generated multiple-sequence alignments with more conserved lengths and qualities than the Pileup program for the different parameters tested (Fig. 1c and d).
The Clustal W program generated the most conserved and consistent alignment with GOP equal to 100 and GEP equal to 100 and a delay divergent option of 90%. The alignment had an overall similarity of 54.78% and was 453 nucleotides long. Instead, GOP equal to 5 and GEP equal to 0.3 or 1 generated the best alignment with the Pileup program, with an overall similarity of 54.61% and a length of 452 nucleotides. Other combinations of parameter values generated very different results.
(ii) Homology index as method of EV identification. The results obtained by use of the homology index are shown in Table 1. The homology test with the Clustal W program identified every serotype but an EV13 strain (homology indices with prototype strains EV13 Del Carmen and ENV69 Tolucal, 74.2 and 73.5%, respectively) and the untypeable strains with values higher than 75%. The PA88-8412 sequence presented the same problem encountered with the pairwise sequence analysis, showing a second-highest score higher than 75%.
The homology test with the Pileup program also identified every serotype but an EV13 strain, the untypeable strains, and field strain T185. Strain T185, typed as EV29, and was correctly identified by the other methods studied (Table 1), while the score obtained with the Pileup program with its respective reference strain was 30.6%. Again, the PA88-8412 sequence resulted in a second-highest score higher than 75% (the score for EV8 Bryson was 77.5%).
(iii) Analysis of test results over multiple reference strain groups. Neither the homology index obtained with the Clustal W program nor the homology index obtained with the Pileup program was affected by the multiple reference strain group factor since all the field strains tested correlated with their corresponding reference strains.
Phylogenetic methods. All the phylogenetic methods generated the same tree topology with different bootstrap values, and the one with higher bootstrap values was obtained with the KITSCH program (Fig. 2). Bootstrap values for each serotype cluster are shown in Table 1. Four separate groups supported by high bootstrap values clustered similarly to other previously published studies (genogroups A, B, C, and D) (15, 20–22). The method identified all serotypes except an EV13 isolate (VA86-6776). This isolate was correctly grouped with its reference strain with a not significant bootstrap value (42%), while it was grouped with ENV69 and EV13 reference strains with higher values. An untypeable EV serotype that could not be identified by other methods of analysis (19) was correctly grouped into genogroup B of the EV genus with a high bootstrap value.
|
|
DISCUSSION |
|---|
|
TopAbstractIntroductionMaterials and MethodsResultsDiscussionReferences |
|---|
The complexity of this genus makes the correct identification of EVs difficult, and problems associated with the seroneutralization method or analysis of untypeable EV strains appear frequently. Nevertheless, although identification of EVs through analysis of the VP1 genome region has been recognized as the best choice, it is also necessary to achieve consensus on a common system of sequence analysis.
The alignment of two sequences is critical for all sequence analysis methods. Differences in alignments and accuracies of methods have previously been reported with large sets of protein sequences with several programs (12). However, the biological validities of nucleotide sequence alignments have been studied less. This is because DNA sequences (which have 4 character states) are more difficult to align than amino acid sequences (which have 20 character states). It was also demonstrated for rRNA data and protein-encoding sequences that alignment accuracy depends on the alignment program and the parameters chosen (14). As was expected, the more that insertion or deletion events are incorporated, the greater the differences among sequences alignments are.
Like other steps in phylogenetic reconstruction, automated sequence alignment requires selection of the most appropriate parameters. Therefore, published studies should describe the parameter values used to create the alignments if they are to be repeatable. In addition, each method of analysis should provide an assessment of quality.
Pairwise nucleotide sequence alignment methods were very different. Although both programs use the same algorithm, the results obtained with the DNASTAR program were very poor.
Pairwise sequence alignment methods provide a simple way to assess the alignment quality by the Monte Carlo approach. The software developer indicates that the algorithm used in the MegAlign program is a slight modification (11) of the algorithm originally designed by Needleman and Wunsch (16) for proteins. This allows the user to align two nucleotide sequences, first by an approach described by Martinez (11), which identifies regions with perfect matches, and then by the method of Needleman and Wunsch (16), which optimizes the fit between perfect matches. The MegAlign program does not provide a quality assessment option, and the results obtained were different from those obtained with the Gap program. An example of this is for isolate R100, which was correctly identified by the Gap program, while the MegAlign program showed a very low similarity index (46.2%) in tests with the same reference strain. In addition, the M-NW similarity test score resulted in the worst performance, and only 32.7% of the strains fulfilled the criteria for identification. Therefore, the MegAlign program is useless for the identification of EVs by use of this genome region.
Instead, the GCG package provides an option that assesses the quality of the alignment obtained with the Gap program. The Monte Carlo simulation determined that the best pair of parameter values was GOP equal to 50 and GEP equal to 3. Other conditions resulted in scores for quality controls lower than the proposed 5 SDs for several pairs of alignments.
By use of these parameters, the highest scores were higher than 75% for 99 of 104 (94.6%) clinical strains. However, 46 clinical isolates could not be identified since they did not fulfill the second-highest-score criterion. Added to that isolate PA88-8412 had a score higher than 75% for two different serotypes, including the correct one (EV4). Scores higher than 70% slightly decreased as the GOP increased (data not shown), showing that these pairwise sequence alignments were relating divergent pairs of sequences, as suggested previously (17). Although the parameters had to be varied to increase the penalty to the similarity score for divergent sequences, the distances between similarity scores were not enough to increase the number of identifiable isolates.
In summary, this pairwise sequence alignment method did not allow the identification of 52 of 104 (50%) strains mainly due to problems related to the proposed criteria. It seems that the criterion of a second-highest score lower than 70% established for identification (19) is very stringent, since all the strains that did not fulfill this criterion had highest scores greater than 75% and the correct serotype.
An even more disappointing aspect was that the pattern of gaps in multiple-sequence alignments was often inconsistent with that in pairwise sequence alignments; thus, the optimal alignment between the two closest sequences was often altered in the presence of a third or a fourth sequence.
Analysis of deduced protein sequences with the Gap program improved the results. Unfortunately, the four untypeable isolates also could not be identified, but the EV13 field strain was correctly identified. The multiple-reference-strain effect was not observed, and the quality analysis had higher values. In this case, the second-highest-score criterion should be also discarded since 25 field strains would have been correctly identified if this criterion had not been considered. The analysis of deduced amino acid sequences could constitute a better approach than the analysis of nucleotide sequences.
The greater precision of the Clustal W program and the more consistent performance of the Clustal W program over that of the Pileup program is due to more sophisticated features such as neighbor joining that guide the alignments. Based on the same reasons considered for pairwise sequence alignment, two different methods were used to evaluate the accuracies of the alignments obtained.
Since length is one of the parameters that allows the user to test the quality of the alignment at first sight, the deduced amino acid sequences were evaluated by considering that if the overall length of a protein alignment is known, the nucleotide sequence alignment length could be calculated (overall amino acid length times 3). The results were different for the two methods. The alignments obtained with the Clustal W program were more consistent and conserved, and better methods to improve quality with the delay divergent sequence option were provided. Moreover, both methods of quality control, SPS and CS, were more consistent for the Clustal W program than for the Pileup program. The alignments of the nucleotide sequences obtained with the Clustal W program were also more reliable than those obtained with the Pileup program. Although the Clustal W program generated the exact 453-bp alignment, the Pileup program was not able to generate an alignment with the predicted overall length. Despite that, the most important features of the Clustal W program are its reproducibility and reliability, given that the resulting alignments were similar under all conditions selected.
The Clustal W program homology index identified all except six clinical strains; while the Pileup program homology index did not identify seven strains. The second-highest-score criterion should not be considered for these tests since their alignments are clearly more reliable than pairwise sequence alignments. In summary, the Clustal W program homology index was able to identify 98 of 104 (94.2%) field strains. Added to that, the Clustal W program is freely available on the Internet and can be performed with a wide variety of computer platforms or with a common personal computer in less than a half hour, while the GCG program is an expensive commercial package that must be run on a mainframe computer system and some Unix programming knowledge is needed to perform the process automatically. This should be considered when proposing a method easily transferable to all reference laboratories.
In any phylogenetic study, alignment of a less conserved region and regions with variable lengths is often problematic, even when sequence divergence is moderate. The alignment process is critical in phylogenetic analyses since inclusion of ambiguously aligned regions could result in erroneous patterns of branching, while removal of such regions could reduce the resolution. To avoid this reduction in the resolution, the terminal regions of the VP1 sequences were included, because it is known that the VP1 region is very divergent and generates a variable-length alignment.
Phylogenetic reconstruction allowed not only the correct clustering of all serotypes with a high bootstrap value (except for the EV13 field strain) but also the prediction and clustering of the new serotypes. This is shown for the cluster of untypeable isolates from clinical samples. Despite the location of the untypeable serotype near the EV20 JV1 reference strain, the cluster is clearly composed of strains of a new serotype since the observed nucleotide sequence distance among strains of the untypeable serotype is lower than 0.20 and the distance between strains of the untypeable serotype and EV20 strains is higher than 0.45 (data not shown).
Identification of the EV13 field strain is a special case since it was clustered with both EV13 and enterovirus type 69 (ENV69) reference strains with high bootstrap values. It was also the only strain that could not be correctly identified by all the methods tested. This could mean that the three strains are of the same serotype, but this hypothesis should be further evaluated with more sequence data from field strains of the EV13 or ENV69 serotype.
In addition, 16 serotypes of EV were not represented in the data. Field isolates of these serotypes are not available in GenBank or in the set of samples sequenced so far. The results should be evaluated again when fields samples of these serotypes are available.
The results of the present study indicate that phylogenetic inference is the best method. Multiple-sequence alignment programs or pairwise sequence alignment of deduced amino acid sequences appears to be the second choice, although adoption of universal criteria for selection of parameters and assessment of the quality of the alignment is needed before multiple-sequence alignment programs can be used. The parameters selected should be also reported in order to repeat and compare the experiments.
Finally, the present study shows that accurate identification of EVs can easily be achieved directly from clinical samples with freely available software.
|
ACKNOWLEDGMENTS |
|---|
The work was supported by institute funds from INEI-ANLIS "Dr. Carlos G. Malbrán" and by Fondo de Investigaciones Sanitarias grant FIS98/0229 from the Spanish Ministry of Health. I.C. is a postdoctoral fellow funded by the Instituto de Salud Carlos III, Becas de Perfeccionamiento.
|
FOOTNOTES |
|---|
|
REFERENCES |
|---|
|
TopAbstractIntroductionMaterials and MethodsResultsDiscussionReferences |
|---|
This article has been cited by other articles:
| |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| Antimicrob. Agents Chemother. | Clin. Microbiol. Rev. |
|---|---|
| Clin. Vaccine Immunol. | ALL ASM JOURNALS |
|---|
