extfieldsofstudy: Medicine Biology

text:
Quantitative Analysis of Severe Acute Respiratory Syndrome (SARS)-associated Coronavirus-infected Cells Using Proteomic Approaches

We present the first proteomic analysis on the cellular response to severe acute respiratory syndrome-associated coronavirus (SARS-CoV) infection. The differential proteomes of Vero E6 cells with and without infection of the SARS-CoV were resolved and quantitated with two-dimensional differential gel electrophoresis followed by ESI-MS/MS identification. Moreover isotope-coded affinity tag technology coupled with two-dimensional LC-MS/MS were also applied to the differential proteins of infected cells. By combining these two complementary strategies, 355 unique proteins were identified and quantitated with 186 of them differentially expressed (at least 1.5-fold quantitative alteration) between infected and uninfected Vero E6 cells. The implication for cellular responses to virus infection was analyzed in depth according to the proteomic results. Thus, the present work provides large scale protein-related information to investigate the mechanism of SARS-CoV infection and pathogenesis.

Quantitative Analysis of Severe Acute Respiratory Syndrome (SARS)-associated Coronavirus-infected Cells Using Proteomic Approaches
We present the first proteomic analysis on the cellular response to severe acute respiratory syndrome-associated coronavirus (SARS-CoV) infection. The differential proteomes of Vero E6 cells with and without infection of the SARS-CoV were resolved and quantitated with twodimensional differential gel electrophoresis followed by ESI-MS/MS identification. Moreover isotope-coded affinity tag technology coupled with two-dimensional LC-MS/MS were also applied to the differential proteins of infected cells. By combining these two complementary strategies, 355 unique proteins were identified and quantitated with 186 of them differentially expressed (at least 1.5-fold quantitative alteration) between infected and uninfected Vero E6 cells. The implication for cellular responses to virus infection was analyzed in depth according to the proteomic results. Thus, the present work provides large scale protein-related information to investigate the mechanism of SARS-CoV infection and pathogenesis.

Molecular & Cellular Proteomics 4:902-913, 2005.
A new type of coronavirus was reported as the causal agent of severe acute respiratory syndrome (SARS) 1 in April 2003, and the genome of the SARS-CoV was sequenced by several groups (1)(2)(3). The properties of the SARS-CoV genome was analyzed in depth by bioinformatic tools (4,5). In addition, several important works on the proteins of SARS-CoV have been reported recently, including the identification of SARS-CoV 3C-like protease structure (6) and the identification of angiotensin-converting enzyme 2 as a functional receptor for the spike protein (7). In our recent work, we identified all of the predicted SARS-CoV structural proteins, nucleocapsid (N), membrane (M), spike (S), and envelope (E), using proteomic approaches and found a novel protein, SARS-CoV 3a (8,9).
To uncover the mechanisms of cellular responses to the virus infection and identify potential drug targets of antiviral treatment, it is very useful to study the molecular profiling of virus-infected cells with high throughput and quantitative approaches. Analysis of gene expression profiles during viral infection is one of the powerful approaches to probe potential cellular genes involved in viral infection and pathogenesis (10). The recent development of proteomic analytic technology such as differential gel electrophoresis (DIGE) (11,12) and ICAT (13,14) also provides new tools for such studies.
As a method based on two-dimensional (2D) electrophoresis, DIGE allows two or three independent samples labeled with different fluorescent dyes such as cyanine-2 (Cy2), cyanine-3 (Cy3), and cyanine-5 (Cy5) to be run in one gel simultaneously and viewed individually using the different fluorescent properties of Cy2, Cy3, and Cy5, circumventing some of the reproducibility problems associated with 2D electrophoresis and providing more accurate quantitative information compared with other staining methods such as silver staining with the dynamic range over 3-4 orders of magnitude (11,12). More recently, the combination of stable ICAT, LC, and MS/MS has emerged as an alternative gel-free quantitative proteomic technology (13,14). In ICAT analysis, two pools of proteins are labeled respectively with isotopically light and heavy ICAT regents, which are chemically identical and therefore serve as a good internal standard for accurate quantification. Although LC-MS provides quantitative information based upon the relative abundances of the heavy and light peptides, LC-MS/MS provides qualitative information based upon the peptide molecular mass and amino acid sequence information. These two technologies have been proved to be complementary for a comprehensive comparative proteomic analysis (15,16). ICAT analysis showed a clear bias for proteins with high molecular mass, whereas the 2D electrophoresis or DIGE method could separate proteins in certain low molecular mass ranges and also identified cysteine-free proteins that were transparent to the ICAT analysis. Moreover ICAT analysis quantifies the sum of the protein species of one gene product, whereas the 2D electrophoresis or DIGE method quantifies at the level of resolved protein species, including post-translationally modified and processed polypeptides.
In the present work, systematic analyses of the proteome of SARS-CoV-infected cells were performed using 2D-DIGE followed by electrospray mass spectrometry identification and ICAT technology coupled with 2D-LC-MS/MS. Using these two complementary methods, the differentially expressed proteome profiles between Vero E6 cells with and without virus infection were created, and a total of 355 proteins or protein spots were identified and quantitated. Further analysis of these data provides the clues for understanding the infection and pathogenesis of the SARS-CoV and the virus-host interactions.

EXPERIMENTAL PROCEDURES
Materials-Analytical reagent grade chemicals were used throughout unless otherwise stated. Water was prepared using a Milli-Q water purification system (Millipore, Bedford, MA). Chemicals used for gel electrophoresis were purchased from Bio-Rad. Cy2, Cy3, Cy5, Precast IPG dry strips, pH 3-10 non-linear, and ECL PLUS Western blotting detection regents were purchased from Amersham Biosciences. Formic acid and guanidine hydrochloride were obtained from Sigma. HPLC grade ACN was from Fisher. Sequencing grade trypsin was obtained from Promega (Southampton, UK).
For virus infection, Vero E6 cells were treated with the DMEM (2% FBS) containing SARS-CoV virions (BJ-01 isolate, provided by Academy of Military Medical Sciences) for 1 h. The virus medium was removed after the infection, and the infected cells were cultured in DMEM with 2% FBS at 37°C in 5% CO 2 . 24 h postinfection, the cells were washed with PBS twice and then lysed with a lysis buffer containing 40 mM Tris base, 60 mM DTT, 8 M urea, 4% CHAPS. The mixture was centrifuged in a microcentrifuge for 5 min, and the supernatant was collected as the infected cellular sample. All the experiments using the virus were carried on in a biosafety Level 3 laboratory.
In-gel Trypsin Digestion-Changes observed in 2D-DIGE images were aligned with Phastgel Blue R-stained protein patterns in preparative gels. Spots of interest were cut manually. Gel pieces were digested as described by Yu et al. (17).
RP-HPLC-ESI-MS/MS Analysis-RP-HPLC was performed using a surveyor LC system (Thermo Finnigan, San Jose, CA) on a C 18 column (RP, 180 m ϫ 150 mm, BioBasic® C 18 , 5 m, Thermo Hypersil-Keystone). The pump flow rate was split 1:120 to achieve a column flow rate of 1.5 l/min. Mobile phase A was 0.1% formic acid in water, and mobile phase B was 0.1% formic acid in acetonitrile. The tryptic peptide mixtures were eluted using a gradient of 2-98% B over 60 min.
The mass spectral data shown here were acquired on an LTQ linear ion trap mass spectrometer (Thermo Finnigan) equipped with an electrospray interface operated in positive ion mode. The temperature of the heated capillary was set at 170°C. A voltage of 3.4 kV applied to the ESI needle resulted in a distinct signal. Normalized collision energy was 35.0. The mass spectrometer was set so that one full MS scan was followed by three MS/MS scans on the three most intense ions from the MS spectrum with the following Dynamic Exclusion TM settings: repeat count, 2; repeat duration, 0.5 min; exclusion duration, 2.0 min.
ICAT Analysis-ICAT analysis was performed using the Cleavable ICAT TM reagent kit (Applied Biosystems, Foster City, CA) according to the manufacturer's guidelines with some modifications. For ICAT analysis, the cellular samples were precipitated and resolubilized in denaturing buffer (6 M guanidine hydrochloride, 100 mM Tris-Cl, pH 8.3). 100 g of the E6 or E6-V protein sample in 80 l of denaturing buffer were reduced at 37°C for 2 h with 5 mM tributylphosphine (Bio-Rad) and alkylated at 37°C for 2 h in the dark with ICAT-light and ICAT-heavy reagent, respectively. The samples were digested with trypsin at 37°C for 20 h. Then the ICAT-labeled peptides were purified using the kit of ICAT TM Avidin Buffer Pack and Avidin Affinity Cartridge (Applied Biosystems) according to the manufacturer's guidelines.
2D LC-MS/MS-Orthogonal 2D LC-MS/MS was performed using a ProteomeX work station (Thermo Finnigan). The system was fitted with a strong cation exchange column (320-m inner diameter ϫ 100 mm, DEV SCX, Thermo Hypersil-Keystone) and two C 18 RP columns (180 m ϫ 100 mm, BioBasic C 18 , 5 m, Thermo Hypersil-Keystone). The salt steps used were 0, 25, 50, 75, 100, 150, 200, 400, and 800 mM NH 4 Cl synchronized with nine 140-min RP gradients. RP solvents were 0.1% formic acid in either water (A) or acetonitrile (B). The setting of the LCQ Deca Xplus ion trap mass spectrometer is as follows. One full MS scan was followed by three MS/MS scans on the three most intense ions from the MS spectrum according to the following Dynamic Exclusion settings: repeat count, 1; repeat duration, 0.5 min; exclusion duration, 3.0 min.
Data Base Searching-The acquired MS/MS spectra were automatically searched against the nonredundant human protein data base (NCBI (www.ncbi.nlm.nih.gov), December 4, 2003 release) using the TurboSEQUEST program in the BioWorks TM 3.1 software suite. For protein spot identification, an accepted SEQUEST result had to have a ⌬Cn score of at least 0.1 (regardless of charge state). Peptides with a ϩ1 charge state were accepted if they were fully tryptic and had a cross correlation (Xcorr) of at least 1.8. Peptides with a ϩ2 charge state were accepted if they had an Xcorr Ͼ2.5. Peptides with a ϩ3 charge state were accepted if they had an Xcorr Ͼ3.7. For ICAT analysis, protein identification and quantification were achieved by using SEQUEST and EXPRESS software tools. Peptides with a ϩ1 charge state were accepted if they were fully tryptic and had an Xcorr of at least 1.5. Peptides with a ϩ2 charge state were accepted if they had an Xcorr Ͼ2.0. Peptides with a ϩ3 charge state were accepted if they had an Xcorr Ͼ2.5. Then the confirmation of protein identification and quantification of the peptides was further analyzed manually as described by Han et al. (14). Protein abundance ratios larger than ϩ1.5 or smaller than Ϫ1.5 (or 0.67) were set as a threshold indicating significant changes (18).
Bioinformatic Annotation-The theoretical pI and molecular mass values of proteins were defined by the program pepstats (www. hgmp.mrc.ac.uk/Software/EMBOSS). The protein function and subcellular location annotation was from the Swiss-Prot and TrEMBL protein data base (us.expasy.org/sprot/).

2D-DIGE Analysis of the SARS-CoV-infected and Uninfected Vero E6 cells-2D
-DIGE as a qualitative and quantitative proteomic approach was performed to determine the differential proteomes of the SARS-CoV-infected and uninfected Vero E6 cells. The lysate from infected cells was labeled with Cy5, whereas the lysate from uninfected cells was labeled with Cy3. After electrophoresis and imaging, the Cy3 and Cy5 images were false colored in red and green, respectively, and two images were overlapped (Fig. 1). The biological FIG. 1. Cy dye image of the SARS-CoV-infected Vero E6 (red, Cy3-labeled) cells and uninfected Vero E6 cells (green, Cy5-labeled). The image shown is of a 100-g sample (50 g each of Cy3-and Cy5-labeled lysates) run on a pH 3-10 non-linear gradient IPG strip and 12.5% polyacrylamide gel. The red or green spots indicate the differentially expressed proteins. The spots circled and marked with numbers have been identified and listed in Supplemental Table I.  Table I). Among those proteins, 17 proteins (21 protein spots) were down-regulated, whereas 13 proteins (15 protein spots) were up-regulated by at least 1.5-fold in SARS-CoV-infected cells (Table I). The identified protein spots with at least a 2.0-fold quantitative alteration in SARS-CoV-infected Vero E6 cells are listed in Table II. The differentially expressed proteins are involved in various functions including cytoskeleton, actin-associated network, metabolic enzymes, and signal transduction (Table II and Supplemental Table I).
ICAT Analysis of the SARS-CoV-infected and Uninfected Vero E6 Cells-Recently a new proteomic approach, ICAT, has been developed to quantitatively analyze the protein differential expression, which is faster than the gel method. In the present work, the ICAT method coupled to 2D-LC-MS/MS was used for quantitative comparison of differential proteome profiles between SARS-CoV-infected cells and uninfected cells. Cysteines were labeled with light ( 12 C) and heavy ( 13 C) ICAT reagent, and labeled peptides were affinity-purified using an avidin column. The resultant peptides were first frac-   a Select list of the peptides with the best assignment scores for protein identification and quantification by ICAT analysis. C/C*, ICAT light/heavy reagent-labeled cysteine-including peptide, respectively. b For heterogeneous nuclear ribonucleoproteins, ribosomal proteins, and translation factors, proteins with at least 1.5-fold alteration are listed. The MS and MS/MS spectra for identification and quantitation of these proteins are shown in Supplemental Figs. 4 -24. tionated with nine salt steps, and each step was further separated using a reverse phase column and analyzed with ESI-MS/MS. The relative quantitation of proteins in infected and uninfected cells was determined from the relative abundance ratio of labeled peptides. By this approach, the relative abundance of 322 proteins was measured (Supplemental Table I).
Individual peptide sequences and their probability scores for identification and quantitation of those proteins are listed in Supplemental Table II. Among these cellular proteins identified from the infected cells, 167 proteins had at least 1.5-fold quantitative alterations with 119 proteins up-regulated and 48 proteins down-regulated (Table I). Some selected proteins with quantitative alterations are listed in Table III. In addition, two SARS proteins, M and S, were observed to be significantly increased in the infected cells (data not shown).
Comparison and Validation of the Quantitation of Differential Proteins-We used 2D-DIGE and ICAT methods to investigate the differentially expressed proteins in Vero E6 cells infected with SARS-CoV. These two methods are both more sensitive and accurate than traditional methods such as silver staining (data not shown). We further used Western blot assay to confirm the differential expression of the proteomes identified by those two methods. The Western blotting results showed that the ratios of four representative proteins, ␣-tubulin, HSP90, OP18 stathmin, and vimentin, between the infected and uninfected cells were in agreement with those obtained from 2D-DIGE or ICAT approaches ( Fig. 2 and Supplemental Figs. 1-3).
It was observed that the overall expression alterations were similar with both quantitative proteomic approaches, although only 15 differential proteins were identified by both 2D-DIGE and ICAT approaches (Table IV). In addition, the DIGE method detected more isoforms of proteins than the ICAT approach did (Table IV).
Subcellular Location and Function Classification of the Differential Proteins-In the present work, a total of 355 unique gene products of SARS-CoV-infected and uninfected cells were identified and quantitated with either the 2D-DIGE or ICAT method of which 186 proteins had at least 1.5-fold quantitative alterations (Table I and Supplemental Table I). Among those 186 differentially expressed proteins, 60 proteins were down-regulated in infected cells, whereas 126 proteins were up-regulated. These proteins were further classified according to their subcellular locations (Fig. 3). Interestingly, the up-regulated proteins in infected cells were mainly located in the nuclei (about 25%, see Fig. 3A), whereas downregulated proteins distributed within the cells (Fig. 3B).
Those 186 proteins with at least 1.5-fold differential expression were further classified according to their functions (Supplemental Table I). Fig. 4 presents the functional categories of these proteins. In these identified SARS-CoV-infected cellular proteins, the enzymes, signal proteins, ribosomal proteins, and heterogeneous nuclear ribonucleoproteins (hnRNPs) consisted of about 45% of the total up-regulated proteins (Fig.   4A). On the other hand, more than half of the down-regulated proteins of the infected cells were involved in the enzymatic reactions, signal transduction, immune responses, and actin networks (Fig. 4B). DISCUSSION From the literature, very few studies have been performed to analyze the interaction between coronavirus and host cells with proteomic approaches. In the present work, we used two quantitative proteomic assays, DIGE and ICAT, to determine the differentially expressed protein profiles of SARS-CoVinfected and uninfected cells. It was noted that only a few proteins were identified by both DIGE and ICAT approaches, suggesting that these two methods are complementary to each other. In addition, the main difference of those two assays is that ICAT only determines the overall expression level, whereas DIGE can detect and quantitate protein isoforms in a gel. Although the DIGE method can detect protein isoforms possibly caused by post-translational modifications, it has limitations in identification of proteins with very high or low molecular weights, extreme acidic/alkaline proteins, and low abundance proteins. In our present work, 45 proteins that have a molecular mass greater than 100 kDa were identified among which only two proteins were detected by the DIGE approach. Moreover all 22 identified proteins with pI Ͼ 10 were contributed by the ICAT method (Supplemental Table I).
Concerning the capacity of the protein identification, the DIGE approach identified 48 unique proteins, whereas the ICAT approach obtained 322 proteins; the ICAT approach especially identified more signal proteins, which usually are low abundance proteins in a cell.
Based on the identified proteins in the present work, we can gain an overall insight into the altered protein expression of the host cell responding to SARS-CoV infection. Among the differentially expressed host proteins, many of them participated in viral RNA replication and translation (Tables II and III and Supplemental Table I). It is known that positive-strand RNA viruses recruit normal components of host cellular RNA-processing or translation machineries for the viral RNA synthesis and protein synthesis (19,20). Our data here indicate that SARS-CoV virus uses a similar strategy when infecting the host cells. hnRNPs are described as a major group of nuclear RNAbinding proteins that function in transcription, RNA processing, mRNA translation, and turnover (21,22). In the present work, a total of eight kinds of hnRNP factors were identified to be up-regulated significantly in SARA-CoV-infected cells ( Table III  requires the function of hnRNPs. Among these identified hnRNPs, hnRNP A1, hnRNP K, and poly(rC)-binding protein have been reported previously to participate in positive-strand virus genome replication. Moreover, the rest of the hnRNPs were described here for the first time to be involved in coronavirus infection (Table III).
HnRNP A1 has been extensively studied for its role in viral RNA replication. Some reports showed that hnRNP A1 could bind the RNA of a mouse coronavirus (mouse hepatitis virus (MHV)) at the 3Ј-end of both plus and minus strands and modulate MHV RNA synthesis (21)(22)(23). However, a recent work argued that hnRNP A1 might not be necessary for MHV viral genome replication or transcription in vivo because the absence of hnRNP A1 in infected cells had no effect on the production of infectious MHV (24). In addition, another experiment revealed that a mouse erythroleukemia cell line, CB3, did not express hnRNP A1 but still supported MHV replication, whereas hnRNP A2/B1, hnRNP A/B, and hnRNP A3 could replace hnRNP A1 in cellular functions and viral infection (25). In the present study, hnRNP A1 has been identified and quantitated by ICAT analysis based on doubly charged and both ICAT light reagent-and heavy reagent-labeled peptide K.YHTVNGHNCEVR.K, a typical trypsin-digested peptide fragmentized at the lysine or arginine carboxyl end, which is marked with the period (Supplemental Fig. 4) and shows up-regulated expression (E6:E6-V, 1:1.90) in SARS-CoV-infected cells (Table III). The present work provides a new explanation for such controversial results on hnRNP A1, i.e. several different hnRNP factors in the infected cells may form a functional hnRNP complex participating in viral RNA metabolism in which one hnRNP factor can be substituted by another without disruption of the function of the hnRNP complex.
The viral genomic RNA of all positive-strand RNA viruses need to be translated by recruited host factors (26,27). A recent report showed that the decrease of 60 S ribosome protein levels reduced a positive-strand virus (Brome mosaic virus)-directed expression in yeast cells (27). The quantitatively proteomic approaches used here revealed that about nine kinds of ribosomal proteins, including components of both 40 and 60 S ribosomal subunits, were up-regulated significantly (Table III and Supplemental Figs. [12][13][14][15][16][17][18][19][20], suggesting that overall up-regulation of the ribosomal protein expression is required for positive-strand RNA virus propagation in host cells.
Translation factors have been well documented to participate in the processes of the virus RNA and protein synthesis (26,28). For example, the elongation factors EF-Tu and EF-Ts were found to bind tightly to the viral RNA-dependent RNA polymerase (29). We also identified the up-regulation of EF-Tu expression in SARS-CoV-infected cells (Table III and Supplemental Fig. 21). Some studies showed that eukaryotic translation initiation factor eIF-4 protein complex is involved in viral protein synthesis (26). On the other hand, the present studies revealed that the expression of eIF-1, eIF-5, and eIF-5A2 was increased in the infected cells, suggesting that these factors play a role in the process of SARS-CoV viral translation (Table  III and Supplemental Figs. [22][23][24].
In summary, we present the first quantitative proteomic work on the cellular responses to SARS-CoV infection, establishing so far the most comprehensive differential proteomic index for SARS-CoV-infected cells. The identified differential profile derived from the infected cellular proteins gives the implications for the infectivity and pathogenesis of SARS-CoV and provides a valuable resource for diagnosis, drug development, and clinical treatment for SARS. * This work was supported by the National Key Basic Research Program. The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked "advertisement" in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.


extfieldsofstudy:
Computer Science
Medicine
Biology

s2fieldsofstudy:
Biology
Computer Science

text:
Clustering protein sequences with a novel metric transformed from sequence similarity scores and sequence alignments with neural networks

Background The sequencing of the human genome has enabled us to access a comprehensive list of genes (both experimental and predicted) for further analysis. While a majority of the approximately 30000 known and predicted human coding genes are characterized and have been assigned at least one function, there remains a fair number of genes (about 12000) for which no annotation has been made. The recent sequencing of other genomes has provided us with a huge amount of auxiliary sequence data which could help in the characterization of the human genes. Clustering these sequences into families is one of the first steps to perform comparative studies across several genomes. Results Here we report a novel clustering algorithm (CLUGEN) that has been used to cluster sequences of experimentally verified and predicted proteins from all sequenced genomes using a novel distance metric which is a neural network score between a pair of protein sequences. This distance metric is based on the pairwise sequence similarity score and the similarity between their domain structures. The distance metric is the probability that a pair of protein sequences are of the same Interpro family/domain, which facilitates the modelling of transitive homology closure to detect remote homologues. The hierarchical average clustering method is applied with the new distance metric. Conclusion Benchmarking studies of our algorithm versus those reported in the literature shows that our algorithm provides clustering results with lower false positive and false negative rates. The clustering algorithm is applied to cluster several eukaryotic genomes and several dozens of prokaryotic genomes.

Background
Clustering of protein sequences from different organisms has been used to identify orthologous and paralogous protein sequences, to find protein sequences unique to an organism, and to derive the phylogenetic profile for a cluster of protein sequences. These are some of the essential components of a comparative genomics study of protein sequences across several genomes.
The methods of clustering protein sequences can be either domain-based or family-based. All the clustering methods start with an all-against-all pairwise protein sequence similarity search. The domain-based clustering methods organize the protein sequence universe into domain clusters where domains are the structural units of proteins, e.g., COG [1], ProDom [2], and Picasso [3] (Figure 1). A multidomain protein may belong to multiple domain clusters.
Clusters of Orthologous Groups (COGs) find triangles of mutually consistent genome-specific best hits from distant organisms without specifying a fixed similarity cutoff, thus accommodating both fast evolving and slow evolving genes. It then merges triangles which share a common edge. Each COG cluster is further analyzed manually to eliminate false positives caused by multidomain proteins so that each COG cluster represents a domain.
ProDom is based on the assumption that short protein sequences are single domain proteins. It first sorts all the protein sequences according to their lengths. It then undergoes a repetitive process: during each iteration, Pro-Dom chooses as the query sequence the current shortest protein sequence or its internal repeat unit if it has internal repeats, searches the whole protein sequence set with PSI-BLAST [4], builds the sequence profile, and masks segments covered by the sequence profile for multidomain The schematic view of family-based clustering Figure 1 The schematic view of family-based clustering. Figure 1 illustrates a typical example of the clustering of three protein families denoted by the three oval outlines. Family I consists of protein sequences 1 and 2. Family II consists of protein sequences 3, 4, and 5. Family III consists of protein sequence 6 and 7. Domain A is common to families 1 and 2 while Domain B is common to families 2 and 3. proteins or removes the single domain proteins completely covered by the sequence profile. The process iterates until there is no sequence left in the protein sequence set.
Picasso merges pairwise sequence alignments from the initial all-against-all pairwise sequence similarity searches into multiple sequence alignments of closer homologs, and later hierarchically merges these multiple sequence alignments into representative sequence profiles of remote homologs by profile-profile comparisons. The representative sequence profiles may contain sequences of different domain structures, but share at least one domain. Picasso then cuts domains within the representative sequence profiles into individual domain clusters based on the concept of overlapping maximal clusters proposed in SYSTERS [5]. Maximal clusters are clusters not fully contained in any other clusters. Two maximal clusters may have not only the overlapping set of neighbour members but also the unique set of neighbour members to these two maximal clusters. Thus these two overlapping maximal clusters must be of different domain structures sharing at least one domain which corresponds to the overlapping set of neighbour members. Then these two overlapping maximal clusters must undergo domaincutting to be split into individual domains corresponding to closed neighbours, where no member has any neighbour outside of the cluster, from multiple alignments. However, since it is still a challenging problem to precisely pinpoint the structure domain border based on primary sequence information [6,7], the performance of the clustering algorithm will be determined by the accuracy of domain demarcations.
SYSTERS uses each protein sequence in the dataset as a seed sequence and applies the single linkage algorithm with a stringent threshold. Thus, each seed sequence has a cluster associated with it. It then merges all the clusters to maximal clusters. The maximal clusters could be either separate maximal clusters corresponding to single domain protein clusters or overlapping maximal clusters representing clusters having multiple domains, but sharing at least one domain.
ProClust uses a different metric to detect whether the aligned two proteins have similar domain structures. The metric value, which scales from 0 to 1, is the ratio of the raw score of the sequence alignments to the raw score of one of those two sequences aligned to itself. Thus the metric value between two sequences is directional. It assumes that the metric is symmetric if two aligned sequences have similar domain structure and non-symmetric otherwise. It then represents each sequence as a vertex and represents the metric value above the threshold as a directional edge in a directed graph. Each strongly connected component corresponds to a cluster [9]. It was later improved by building Profile-HMMs for all clusters having more than 20 sequences and merging two clusters A and B into a cluster corresponding to a SCOP superfamily if the average Evalue from searching all the sequences in the cluster A against the profile-HMM of the cluster B is below the threshold [10].
PROTONET applies the hierarchical clustering of protein sequences based on the their pairwise similarity E-values, but adopts different rules for merging clusters: arithmetic mean, geometric mean, and harmonic mean. However, different families may have different levels of sequence conservation. It is not appropriate to choose one E-value threshold. And at the level of higher E-value, it may merge two clusters of different domain structures, but sharing one domain. However, different families may have different levels of sequence conservation. It is not appropriate to choose one E-value threshold. And at the level of higher E-value, it may merge two clusters of different domain structures, but sharing one domain.
Transitive homology detection methods have been proposed in the Intermediate Sequence Search, ISS [14,15], and [16]. It works by searching the query sequence against the database with a conservative threshold to find the closely homologous sequences and using these homologous sequences as seeds to search the database to find remotely homologous sequences with a less conservative threshold. The method has been shown to be close to the profile based methods and better than a direct pairwise homology search [17]. But, it is a challenge to quantify the indirect, transitive homology as opposed to using the Evalue for quantifying direct pairwise sequence homology.
The Markov cluster (MCL) [13] algorithm has been successfully applied to clustering protein sequences. MCL represents protein sequences as nodes on a graph where similar proteins are connected by edges weighted according to their BLASP E-Value. The MCL algorithm works through a series of iterative random walks across the graph and inflations of the edge weights that gradually strengthens the connections between very similar nodes and weakens the connections between less similar nodes. MCL makes no explicit use of protein domain architecture but does leverage transitive homologies in the random walk phase of the algorithm.
Compared to the hierarchical clustering family based clustering method, e.g., PROTONET [11], our method can take advantage of the transitive homologue closure by the third intermediate sequence to detect remote homologues at the superfamily level. Compared to single linkage based methods, e.g., CluSTr [8], our method avoids the problem of merging two unrelated multi-domain cluster sharing a common domain. Compared to the iterative clustering method, e.g., SYSTERS [5], our method generates clusters where each sequence belongs to only one cluster.

Benchmarking
In order to test the performance of CLUGEN, we selected all Swissprot [18] sequences with an InterPro [19] annotation, which resulted in 41480 sequences from 1598 InterPro families. The criteria we used to select sequences are that more than one member database from the Inter-Pro annotation have the same superfamily or domain assignment and that the aligned region of the Swissprot sequence with respect to either Profile or hidden Markov model is longer than 30 amino acids. The benchmarking dataset is available on request. Figure 2 shows the InterPro superfamily/domain size distribution in the benchmarking dataset. There are 102 singleton families, that is families that consist of only one sequence. The largest family is IPR000276, the Rhodopsin-like G-Protein Coupled Receptor (GPCR) family which comprised of 1058 protein sequences.

Performance measure
We measure CLUGEN's performance by sensitivity, specificity, and goodness. A protein sequence is a false positive (FP) if it is misclassified to a certain InterPro superfamily/ domain and a true positive (TP) otherwise. A protein sequence in a certain InterPro superfamily/domain is a false negative (FN) if it is not classified to that InterPro superfamily/domain ( Figure 3). Let N fp and N tp denote the number of false positives and the number of true positives with respect to a cluster. Let N fn denote the number of false negatives with respect to a InterPro superfamily/ domain.
shows the distribution of InterPro family size Figure 2 shows the distribution of InterPro family size. Figure 2 shows the distribution of the InterPro families used in the benchmarking dataset based upon the number of members in each family. There are 102 singleton InterPro families, and the largest InterPro family in the benchmarking dataset is Rhodopsin-like GPCR superfamily which has 1058 protein sequences in the benchmarking dataset. Goodness: The goodness of an InterPro superfamily/ domain is a measure of how well a cluster corresponds to the mapped InterPro superfamily/domain and is defined below (Equation 1) where N denotes the number of generated clusters associated with that InterPro superfamily/ domain. The Area Under the ROC Curve (AUC) has been shown to be a better evaluation measure than accuracy within the context of binary classification, where the negative dataset is clearly defined. However, we cluster protein sequences into 1598 interpro families simultaneously. As a result, using the AUC as a measure of performance is not the appropriate metric here. Instead, we adopt as the "goodness " the set relative measure as defined in [12]. In order to decrease the goodness value when a large number of clusters is associated with an InterPro superfamily/domain, a penalty of (N-1) is applied in the numerator of the equation.

Equation 1:
Definition of various clustering parameters Figure 3 Definition of various clustering parameters. Figure 3 illustrates the mapping of three generated clusters denoted by oval outlines differentiated by different colors into a InterPro family denoted by a rectangle. The cluster can be mapped to an Inter-Pro family only if more than 50% cluster members belong to that InterPro family; and is declared as a orphan cluster otherwise. Protein sequences outside the rectangle are false positives. Protein sequences within both the oval outline and the rectangle are true positives. Protein sequences wholly within the grey rectangle are false negatives.

Overall performance
We evaluated CLUGEN at several threshold values. Table  1 lists shows the specificity, sensitivity, and goodness as well as the number of generated clusters and the number of orphan clusters as a function of the different threshold values respectively. A cluster can be mapped to an InterPro family only if more than 50% of the cluster members belong to that InterPro family; and is declared as an orphan cluster otherwise. At one extreme of the spectrum, each cluster is a singleton cluster consisting of only one protein sequence. Thus both specificity and sensitivity are 100%. But the goodness value is very low, the reciprocal of the size of the family. Clustering based on more stringent threshold values, e.g. 0.5, generates a larger number of smaller clusters causing a smaller number of false positives, also resulting in a low goodness value. As the threshold values become less stringent, small clusters merge at different levels into larger clusters, therefore a smaller number of larger clusters are generated. At the threshold of 0.2, there are 1706 clusters resulting in a specificity of 97.1%, sensitivity 98.6%, goodness value of 78.2%, and the number of orphan clusters is 201. As can be seen from table 1, the threshold value is a compromise of sensitivity, specificity, goodness and the number of orphan clusters. Ideally, we would like the clustering results to produce as many clusters as there should be and as few orphan clusters as possible.
For a basis of comparison we also applied the MCL [13] algorithm to the same test dataset with various inflation values. Results are depicted in Figures 4 and 5. At higher specificities, the sensitivity of both methods increases. This is expected because higher specificities are achieved via stricter thresholds that result in more clusters overall and fewer large clusters. In the extreme case one could place each test sequence in its own cluster of size 1 and achieve 100% sensitivity and 100% specificity but with a low goodness score. This trade-off between sensitivity, specificity, and goodness is clearly evident in Figure 4; as specificity increases, sensitivity increases whereas goodness decreases.
In Figure 5 we see additional tradeoffs between specificity and overall performance. As specificity increases the number of orphan clusters decreases. This improvement in performance comes with an increase in the total Specificity, sensitivity, and goodness on the benchmarking dataset Figure 4 Specificity, sensitivity, and goodness on the benchmarking dataset. Sensitivity and specificity for CLUGEN and MCL at various specificities. At higher specificities, the sensitivity of both methods increases, whereas the goodness of both methods decreases. This is expected because higher specificities are achieved via stricter parameter thresholds that more clusters overall and fewer large clusters. Performance for both methods is comparable in this range with CLU-GEN performing better at lower specificities and MCL performing better at higher specificities.
number of clusters. Once again the extreme case of one sequence per cluster guarantees no orphan clusters at the cost of many non-informative clusters. Ideally one wishes to strike a balance reducing the number of orphan clusters while not drastically increasing the total number of clusters.
The overall performances of MCL and CLUGEN are fairly similar, with CLUGEN demonstrating a clear advantage at specificities below 0.98. CLUGEN's sensitivity and goodness are better at specificities below 0.98, whereas MCL's goodness is slightly better at specificities greater than 0.98. The number of total clusters and orphan clusters generated by both methods are comparable at specificities below 0.98. CLUGEN tends toward fewer orphan clusters at the cost of more total clusters at higher specificities.

Analysis of some CLUGEN generated clusters
In this section, we will give examples of some successfully generated clusters with one-to-one correspondence to spe- We also conducted a detailed analysis of clusters that had false negatives/false positives in order to understand the areas in which the clustering algorithm could be further improved. The following is a description of errors encountered in clustering algorithms with specific reference to the data from our method.

Errors from low-complexity and coiled-coil regions
The first type of error is due to the presence of low complexity sequences with repetitive sequence patterns or sequences with coiled-coil structures, since we mask the low complexity regions and coiled-coil regions before the all-against-all pairwise similarity searches. As an example, the InterPro family IPR000533, Tropomyosins, which regulate muscle contraction, are alpha-helical proteins that form a coiled-coil. There are 25 tropomyosin sequences in the benchmarking dataset, among which 24 tropomyosin sequences are false negative sequences and appear in the following cluster along with members of IPR002699 ATP synthase subunit D.

Errors from short sequences or from an abundance of certain amino acid type in the sequences
Short sequences with less than 70 amino acids could also cause false positives in the clustering results. Cluster 1259 which is mapped to InterPro family 003019, the metallothionein superfamily, consists of 125 short protein sequences with 68 amino acids on average in length among which 34 false positive protein sequences are from InterPro family IPR001762, Disintegrin, and 34 false positive protein sequences are from InterPro family IPR000877, Bowman-Birk serine protease inhibitor. The reason these families cluster together is that The number of generated clusters and orphan clusters on the benchmarking dataset Figure 5 The number of generated clusters and orphan clusters on the benchmarking dataset. Total clusters and orphan clusters for clugen and MCL at various specificities. With stricter parameter thresholds, overall specificity and the total number of clusters increases for both methods. The larger number of small clusters at higher specificities leads to a reduction in the number of orphan clusters in both methods.
metallothioneins are small proteins with high content of cysteine residues, while disintegrins and Bowman-Birk serine protease inhibitors are also short cysteine-rich protein sequences. This suggests that a more stringent threshold should be applied to cluster short protein sequences which are rich in a particular amino acid.

Similar domain structures in different superfamilies
Sequences that belong to different superfamilies but share similar domain structures may also cluster incorrectly in some cases. For example, 1039 out of 1058 sequences correctly cluster into the IPR000276 Rhodopsin-like G-protein coupled receptors family. But 17 out of 19 sequences from the IPR000276 Rhodopsin-like G-protein coupled receptors which have Cysteine-rich N-terminal regions are mistakenly clustered into the InterPro 000372 which is annotated as a cysteine-rich flanking region N-terminal. Similarly, members of the IPR001878 CCHC Zinc finger domain have been incorrectly clustered into the cluster 1008 which is mapped to the InterPro family 000981 Neurohypophysical hormone because they share the two cysteine residues and other surrounding weak motifs.

Conclusion
This paper describes a novel method for the clustering of protein sequences based on a new metric derived from prediction using neural networks and further utilizing the metric to model the transitive sequence homologue to detect the remote homologue. Good performance with respect to the InterPro protein sequence database has been achieved on the benchmarking dataset.
Clustering of sequences has many applications in target discovery and gene functionalization. One can identify in silico, antimicrobial drug targets by examining clusters without any eukaryotic member sequence in it. These proteins could be potentially selective targets for infectious diseases due to the absence of any appreciable homolog in the human proteome. Such computationally derived targets need to be further validated using experimental data derived from gene expression profiling and proteomics experiments [20]. Another application is to predict the function for prokaryotic proteins of unknown function by phylogenetic profiling [21], where the phylogenetic profile for a cluster is a vector of binary values, with 1 meaning the presence of a genome in that cluster and 0 otherwise. The assumption here is that genes with the same phylogenetic profile could have the same function.

Feature extraction
After we mask the low complexity regions and the coiledcoil regions and carry out the all-against-all pairwise sequence similarity searches, we extract four sets of features to represent the homology between any given pair of sequences. The first two sets of input features detect the homology of two aligned sequences, the last two sets of input features test whether two aligned sequences have similar domain structures. We use neural networks to map input features to a new metric, a probability value which scales from 0 to 1 and is interpreted as the likelihood that two sequences are of the same homologous superfamily.
The first input feature is the log scale of the pairwise Evalue.
The raw score, from which the E-value of the two aligned sequences is derived, is calculated by summing up the log score of each position in the alignment, which assumes that each position is independent of the other. However, in practice, it has been shown that two consecutive positions in the alignments are quite correlated [22]. To model the correlation between two consecutive positions in the alignment, we adopt the concept of the 2-gram encoding method [23]. Ideally, hydrophobic regions of one sequence should align with the hydrophobic regions of the other sequence, and hydrophilic regions should align with each other as well. Each position in the alignment could fall into one of four categories: residue identity denoted by A 1 , hydrophobic similarity denoted by A 2 , hydrophilic similarity denoted by A 3 , and everything else denoted by A 4 . Let Len a denote the total length of the alignment and Occur i,j denote the number of occurrences of i and j, where i is immediately followed by j, with i and j denoting any one of A 1 , A 2 , A 3 , or A 4 respectively. Let freq i,j denote frequency of i,j, which is equal to Occur i,j / (Len a -1). Thus the second set of input features consists of freq i,j values of the alignment positions, which consists of 16 input feature values for a pair of aligned sequences A schematic view of a pairwise alignment Figure 6 A schematic view of a pairwise alignment. Figure 6 shows a pairwise alignment between two aligned sequences. The aligned regions of the two sequences are highlighted. Their boundaries are pinpointed by the arrows.  since each of the two consecutive positions could be one of four possible values.
Each of the two aligned sequences is separated into three segments, the unaligned N-terminal region, the aligned region, and the unaligned C-terminal region, by the beginning position of the alignment, denoted by begin i , and the end position of the alignment, denoted by end i with respect to sequence i ( Figure 6). Let Len i denote the length of sequence i . Then lengths of three segments of sequence i are begin i -1, end i -begin i +1, and Len i -end i , respectively. If we normalize the length of each of three segments within an aligned sequence i by dividing the length of each segment by Len i , we get a vector of three values, Seg i1 = (begin i -1)/Len i , Seg i2 = (end i -begin i +1)/Len i , and Seg i3 = (Len i -end i )/Len i . Intuitively, if the two aligned sequences have similar domain structures, the alignment will divide the two aligned sequences i and j in a similar proportion, and the linear correlation coefficient, LCC 1 defined by Equation 2, between these two vectors tend to be close to 1. So the third set of input features include LCC 1 , Seg i1 , Seg i2 , Seg i3 , Seg j1 , Seg j2 , and Seg j3 .

Equation 2:
The fourth and final input feature is to measure the overlap between two neighbor sets of aligned sequences, where the neighbor set, Set i of sequence i , is defined as the set of protein sequences that sequence i hits when sequence i is used as the query sequence. One straightforward method to measure the overlap is to use the ratio of the cardinality of the intersection of two neighbor sets to the cardinality of the union of two neighbor sets.
Here we propose another method to measure the overlap. Specifically, if we represent the neighbor set of sequence i as Vector i , the value of the k th element of Vector i is the log Evalue, Log ik between sequence i and its k th neighbor in Set i . However, Vector i and Vector j for two aligned sequences, sequence i and sequence j , may be of different dimensions since the cardinalities of Set i and Set j may be different. We make these two vectors have the same dimension by adding the log E-value threshold to Vector i whenever the sequence i has no corresponding neighbor in the neighbor set, Set j of the other aligned sequence j. Thus the last input feature is the linear correlation coefficient LCC 2 between Vector i and Vector j defined by Equation 3. Intuitively, the more similar the domain structure two aligned sequences have, the more similar neighbor sets they will have, and the closer to 1 the linear correlation coefficient will be.

Equation 3:
where N is the dimension of Vector i and Vector j and Log i and Log j are the mean values and are defined by , and .
To summarize, the first input feature is the log scale of the pairwise E-value. The second input feature consists of 16 feature values, which are frequency values of the alignment positions and the third input feature includes 7 values, which relate to the details of the alignments. And finally, the fourth input feature includes 1 feature value which measures the overlap of the neighbour sets to make a total of 25 input features which will be used to train the neural network as described below.

Neural networks
After we represent the sequence homology between a pair of sequences by the set of 25 input features, we train the neural network using the training data. Each homologous pair of sequences is labeled as 1 if they belong to the same InterPro superfamily or the same domain if they are single domain proteins, and 0 otherwise. We selected as large a number of sequences as possible to train the neural networks to avoid overfitting the data. In all, we selected 27844 homologous sequence pairs as the positive training set and 29999 non-homologous sequence pairs as the negative training set.
The neural network we use is a fully connected feed-forward back propagation neural network and has one hidden layer with sigmoid activation functions (see Figure 7). The output layer of the neural network has one output unit. The output value is bounded between 0 and 1. The network is trained with the scaled conjugate gradient algorithm [24] implemented in MATLAB [25].
Given the large amount of training data relative to the number of the weights in the neural network, the neural network is unlikely to overfit. It may however underfit if there are not a sufficient number of weights. If the training data are smaller relative to the number of weights in the neural network, measures should be taken to avoid the overfitting problem and the cross-validation method should be used to choose the best model. Clearly in this study, such is not the case.
We used a split-sample approach in which the validation set is not used during training, but is used to select the best model. After the neural network is trained, it is validated on the validation dataset containing 250597 homologous pairs of sequences and 30000 non-homologous pairs. Different numbers of hidden layer nodes have been tested. The ultimate goal is not to determine if any two proteins sequences belong to the same Interpro family, but to cluster all sequences in Interpro families as accurately as possible. We selected the model with the smallest number of weights and smallest error on the validation set. Thus, we chose 4 hidden layer nodes such that the neural network has the least number of hidden units and the best performance on the validation dataset with a specificity of 94.18% and a sensitivity of 91.81%.

Modeling the transitive homology
The neural network is then used to calculate the metric value for any pair of protein sequences that hit to each other below the E-value that was used as a cutoff. The metric value, P (A,B), for protein sequences A and B is interpreted as the likelihood that these two protein sequences belong to the same InterPro superfamily or have the same single domain. The value P(A,B) is replaced by P(A,C)P (C,B) if there exists a sequence C such that P(A,C)P (C,B) is larger than the current value of P (A,B). This transformation takes advantage of the transitive homology of sequences A and B through the intermediate sequence C, assuming that protein sequences A and C and protein sequences B and C are independently homologous. Figure 8 illustrates the transitive homology between sequence a and sequence b through the third sequence c. The E-values between sequence a and sequence c, sequence c and sequence b, as well as sequence a and sequence b are 0.01, 0.005, 20 respectively. P(a,c), P(c,b), and P(a,b) are 0.8, 0.9, and 0.2 respectively. The homology between sequence a and sequence b cannot be detected with their direct E-value. However, the value of P(a,b) is assigned to 0.72 because of the transitive sequence homology.

Hierarchical average linkage clustering
Once the metric value for every pair of protein sequences is calculated, the hierarchical average linkage clustering method is applied to cluster the protein sequences in the new metric space using the geometric mean as the merging rule.
Hierarchical average linkages uses the Unweighted Pair-Group Average (UPGA) [26] to measure the distance between clusters. Let D i , i = 1,2, ... n. denote the protein sequences contained in Cluster D and let E j , j = 1,2, ..., m denote the protein sequence contained in Cluster E. The geometric mean distance G between Cluster D and Cluster E is defined as Equation 4:

Equation 4:
The hierarchical average linkage clustering works in an iterative process: it begins with each protein sequence as a singleton cluster; during each iteration, it finds two clusters with the lowest metric value, then joins these two clusters into a new cluster, and updates the metric value between this new cluster and all others. This process iterates until the current lowest metric value exceeds the threshold.

Authors' contributions
QM carried out the design and implementation of the method and wrote the manuscript. JDS compared the performances between CLUGEN and MCL. GWC and RC participated in the project. NRN directed and participated in the project and prepared the figures in the manuscript. All authors involved in reviewing and revising the manuscript and approved the final manuscript.


extfieldsofstudy: Biology Medicine

text:
The biology of genomes: sequence gives way to function

A report on the meeting 'The Biology of Genomes', Cold Spring Harbor, USA, 11-15 May 2005.

* To whom correspondence should be addressed. E-mail: lonserr@ninds. nih.gov [ These authors contributed equally to this work.

A B S T R A C T Background
To determine the origin of the neoplastic cell in central nervous system (CNS) hemangioblastomas in von Hippel-Lindau disease (VHL) and its role in tumor formation and distribution, we characterized and differentiated neoplastic cells from hemangioblastomas removed from VHL patients.

Methods and Findings
A total of 31 CNS hemangioblastomas from 25 VHL patients were resected and analyzed. Tumor cells from the hemangioblastomas were characterized, grown, and differentiated into multiple lineages. Resected hemangioblastomas were located in the cerebellum (11 tumors), brainstem (five tumors), and spinal cord (15 tumors). Consistent with an embryologically derived hemangioblast, the neoplastic cells demonstrated coexpression of the mesodermal markers brachyury, Flk-1 (vascular endothelial growth factor-2), and stem cell leukemia (Scl). The neoplastic cells also expressed hematopoietic stem cell antigens and receptors including CD133, CD34, c-kit, Scl, erythropoietin, and erythropoietin receptor. Under specific microenvironments, neoplastic cells (hemangioblasts) were expanded and differentiated into erythrocytic, granulocytic, and endothelial progenitors. Deletion of the wild-type VHL allele in the hematopoietic and endothelial progeny confirmed their neoplastic origin.
Introduction von Hippel-Lindau disease (VHL) is an autosomal dominant neoplasia syndrome resulting from a germline mutation in the VHL tumor suppressor gene on the short arm of Chromosome 3 [1]. VHL has an approximate prevalence of one in 39,000 and over 90% penetrance by 65 years of age [2,3]. Germline mutations in the VHL gene lead to the development of a number of benign or malignant tumors and cysts in multiple organ systems [4]. Affected individuals may develop central nervous system (CNS) tumors, including hemangioblastomas and endolymphatic sac tumors [4][5][6][7]. Visceral manifestations include renal cysts and carcinomas, pheochromocytomas, pancreatic cysts, and neuroendocrine tumors, as well as epididymal and broad ligament cystadenomas [4,[8][9][10].
CNS hemangioblastomas are the most common tumor in VHL, affecting 60%-80% of all patients, and are a frequent cause of neurologic morbidity and mortality [4,7,11]. Hemangioblastomas occur in highly defined regions of the CNS that include the retina, cerebellum, brainstem, and spinal cord [4,7,11]. Histologically, hemangioblastomas are benign, highly vascular tumors composed of neoplastic stromal cells and vessels. Detailed morphologic analyses demonstrate intratumoral blood island formation indicative of extramedullary hematopoiesis [12,13].
Although the anatomic distribution and histologic features of hemangioblastomas in VHL suggest they may be derived from a developmental multipotent cell, the identity and function of the cell of origin in hemangioblastomas remains unknown. To determine the tumor cell of origin, the basis for the unique CNS distribution of hemangioblastomas, and the multipotent capability of hemangioblastomas in VHL, we analyzed these tumors for embryologic and hematopoietic stem cell markers. Based on the characterized features of these tumors, we then derived hematopoietic and endothelial progeny by microenvironment manipulation of resected hemangioblastoma neoplastic cells.

Methods Patients
Twenty-five VHL patients who underwent resection of a CNS hemangioblastoma(s) between 2003 and 2005 at the NIH were included. Patients were screened for the presence of mutations and deletions of the VHL gene. All 25 patients met the diagnostic criteria for VHL [4]. Patients underwent serial craniospinal magnetic resonance imaging (MRI), and the location of all imaged hemangioblastomas was recorded. Research was performed under an institutional review boardapproved protocol (National Institute of Neurological Disorders protocol 03-N-0164) and informed consent was obtained from all patients.

Nested Reverse Transcription Polymerase Chain Reaction
Cells morphologically consistent with a myeloid origin were analyzed for the expression of CD13, a myeloid marker, by nested reverse transcription PCR (RT-PCR) [14]. Briefly, total RNA was extracted and reverse-transcribed and amplified by two rounds of PCR. Product of the secondround PCR was electrophoresed on agarose gel and stained with ethidium bromide.

Loss of Heterozygosity Analysis
Investigated tissues included uncultured tumor cells, cultured tumor cells, cultured nucleated red blood cells, cultured immature white blood cells, and cultured endothelial cells. Peripheral blood-derived mononuclear cells from patients served as a negative control. Specific cell types (identified morphologically and immunohistochemically; as described in text) were individually microdissected from slides for DNA extraction [15], and analyzed by PCR for loss of heterozygosity using the microsatellite markers D3S1038 and D3S1110 (Invitrogen) that flank the VHL gene [13]. The samples were electrophoresed on a polyacrylamide gel for detection of the wild-type VHL allele.

Identification of Tumor Cell of Origin
All tumors were histologically confirmed to be hemangioblastomas by routine staining. Immunofluorescence microscopy demonstrated that the hemangioblastoma neoplastic stromal cells expressed the mesodermal marker brachyury. Brachyury, a founding member of the T-box family of transcription factors, has a conserved role in mesoderm differentiation in vertebrates. Developmental studies indicate that early specification of posterior mesoderm and formation of the notochord is regulated by brachyury [16]. Coexpression of brachyury and Flk-1 (VEGF receptor-2) confirmed that the neoplastic cells were derived from the hemangioblast subset of mesodermal cells ( Figure 2). Flk-1 (VEGF receptor-2) is a type 3 receptor tyrosine kinase that is believed to play a critical role in developmental angiogenesis and hematopoiesis [17]. Expression of Scl, which is essential for both primitive and definitive hematopoiesis, was demonstrated in the tumors (Figure 2). Expression of these embryologic markers was confirmed by Western blot analysis ( Figure 3).

Cellular Markers
To characterize the multipotent capacity of the hemangioblastoma neoplastic cells, we looked for expression of known stem cell markers. Immunofluorescence demonstrated expression of the stem cell marker CD133 ( Figure 2C) and the hematopoietic stem cell marker CD34 in a few scattered cells throughout the hemangioblastomas. Similarly, c-kit and the common myeloid antigen, IL-3 receptor, were expressed in tumor cells in scattered regions throughout the hemangioblastomas. Coexpression of erythropoietin and erythropoietin receptor was demonstrated throughout the tumors ( Figure 2D).

Primary Culture of Neoplastic Cells
At concentrations of 50 U/ml of erythropoietin or more, hemangioblast cells were sustained and expanded in culture ( Figure 4). The identity of the cultured hemangioblasts, which was confirmed by morphology and immunofluorescence, demonstrated coexpression of erythropoietin and erythropoietin receptor. At culture concentrations of erythropoietin at 100 U/ml, two other populations of nucleated cells (nonadherent and adherent) appeared. The nonadherent cells were identified by Giemsa staining as nucleated erythrocyte or granulocyte progenitors in various stages of differentiation ( Figure 4). These hematopoietic progenitor cells expanded over two weeks when the expansion stopped and the absolute number of cells remained stable. Adherent cells were identified as endothelial cells morphologically (flat and elongated) and by expression of endothelial markers Tie-2 and CD31. These cells expanded beginning at two weeks and continued to expand while in culture (at least 90 d).

Origin of Hematopoietic and Endothelial Progeny
To determine their origin, VHL gene deletion analyses were performed on microdissected portions of uncultured tumors, cultured tumor cells, cultured hematopoietic progeny, and cultured endothelial progeny. Deletion of the wildtype VHL allele (loss of heterozygosity) occurred in the uncultured tumor specimens, cultured tumor cells, cultured hematopoietic progeny, and cultured endothelial progeny, confirming their neoplastic origin (Figure 4). Peripheral blood-derived mononuclear cells from patients served as a negative control (Figure 4).

Development and Expansion of Mature Hematopoietic Progeny
To determine if the tumor-derived nonadherent immature hematopoietic progenitors could be expanded into mature erythrocytic or granulocytic progeny, we developed a culture environment designed to expand immature hematopoietic cells ex vivo that exploited the characterized receptor expression (erythropoietin, IL3, and c-kit). In this environment, expansion of mature enucleated erythrocytes and granulocytes occurred by day 3 with maximal proliferation occurring by day 7. In vitro senescence of mature hematologic progeny occurred beyond day 7. The mature erythrocytes displayed characteristics similar to that of normal erythrocytes, including a mean cell volume of 80.5 fl, a mean cell hemoglobin of 29.7 pg, and a mean corpuscular hemoglobin concentration of 37 g/dl. The mature granulocytes displayed characteristic morphology and the presence of CD13 mRNA by RT-PCR.

Expansion of Endothelial Progeny
To expand the tumor-derived immature endothelial progeny, we used a culture media designed to expand immature endothelial cells ex vivo that exploited the characterized receptor expression (erythropoietin, IL3, ckit, and Flk-1 [VEGF receptor-2]). The adherent endothelial cells continued to proliferate while in culture (at least 30 d). These expanded cells were adherent and morphologically similar to endothelial cells (flat and elongated). Consistent with an endothelial lineage, these cells expressed the markers Tie-2 and CD31 ( Figure 5).

Discussion
Identification of the cell of origin in tumors will improve understanding of their pathobiology and treatment. For familial tumor suppressor syndromes, it may also explain the unique tissue distribution of tumor involvement. We found that the neoplastic cells of origin in CNS hemangioblastomas from VHL patients are derived from embryologically arrested mesoderm cells and are hemangioblasts with hematopoietic and vasculogenic (endothelial cell) potential. These cells retain their multipotent differentiation ability and are able to develop into hematopoietic and endothelial progenies. The identification of the cell of origin of this tumor enhances understanding of the biological basis and tumor distribution of this familial tumor suppressor syndrome and may provide therapeutic opportunities.

Cell of Origin for Hemangioblastoma
From analysis of chick embryo development, Sabin (1920) hypothesized the existence of a multipotent embryologic precursor cell capable of forming both blood and vessels [18]. Based on Sabin's hypothesis and the discovery of intratumoral blood and vessel formation in CNS hemangioblasto-  mas, Stein et al. (1960) suggested that an arrested angiomesenchymal cell or embryologic defect was the origin of these tumors [12]. Based on embryonic stem cell lineage studies, Choi et al. (1998) discovered the multipotent embryonic precursor for both hematopoietic and endothelial cells that Sabin had hypothesized, and defined it as a ''hemangioblast'' [17]. Despite these observations, the histogenesis of hemangioblastomas has been unclear. It was only recently that Vortmeyer et al. (1997) identified the stromal cell in hemangioblastomas as the underlying neoplastic cell [19]. These previous findings and recent insights into the molecular basis of embryologic development [16,20,21] and tumor formation now permit investigation into the origin and multipotent capability of the neoplastic cells in hemangioblastomas.
To investigate whether the cell of origin in hemangioblastomas is the embryologically arrested hemangioblast, we characterized the neoplastic stromal cells from CNS heman-gioblastomas derived from VHL patients. Consistent with a cell of mesodermal origin, the neoplastic stromal cells in these hemangioblastomas expressed the protein brachyury [16]. Expression of brachyury, Flk-1 (VEGF receptor 2), and Scl in the neoplastic stromal cells confirm that the tumor cell of origin in hemangioblastomas is an embryologic hemangioblast [20]. Expression of Scl confirms that these are hematopoietic stem cells with self-renewal capabilities and that they are committed to hematopoietic and endothelial lineages [22][23][24][25][26]. These findings also indicate that CNS hemangioblastomas are derived from committed hematopoietic precursors and not a neurogenic progenitor, as has been described with other primary CNS malignancies [27,28].
Because brachyury is expressed only during early mesoderm development [16], its expression in these neoplastic cells suggests that VHL-related hemangioblastomas originate from embryologically arrested mesodermal cells. This finding has several potential implications for hemangioblastoma development in VHL. First, it suggests that the precursor cells (hemangioblasts) for hemangioblastoma formation are present during embryologic development. Second, it also implies that loss of the wild-type allele, which is necessary for hemangioblastoma formation in VHL [29], occurs during embryogenesis. Third, these observations indicate that VHL tumor suppressor gene function affects development. These findings are consistent with previous studies that demonstrate that loss of tumor suppressor gene function can affect normal tissue development [30][31][32]. Recently, Zhu et al. [33] found that mice genetically engineered to lack the neurofibromatosis-1 tumor suppressor gene demonstrated a variety of developmental defects that resulted in increased proliferation of glial progenitors cells, which in some cases led to optic nerve glioma formation (part of the neurofibromatosis-1 neoplasia phenotype). Fourth, an embryologically arrested hemangioblast may explain the presence of fetal hemoglobin within blood islands of resected hemangioblastomas from VHL patients [13,34]. Finally, these data support an emerging concept of tumor formation during development in a variety of different tumor types. Specifically, Trichopoulos et al. [35] has hypothesized a perinatal site-specific increase in stem cells as an underlying mechanism of breast cancer development and Samuelson et al. [36] found that increased head circumference was positively correlated with the incidence of childhood brain cancer, implying that pathogenesis occurs during fetal life.

Development of Hematopoietic and Endothelial Progeny
To determine if the hemangioblasts in hemangioblastomas have the same multipotent developmental capacity (hematopoietic and endothelial) as hemangioblasts during embryologic development, we further characterized these cells and analyzed their capacity for multilineage differentiation. Tumors were characterized by scattered expression of hematopoietic stem cell markers CD133 and CD34, as well as stem cell factor receptor (c-kit) and IL-3 receptor. Based on the uniform expression of erythropoietin and erythropoietin receptor by the hemangioblasts, suggesting a neoplastic erythropoietin-driven autocrine loop [13,37], we initially cultured the neoplastic cells in an erythropoietin-enriched media. Similar to embryologic hemangioblasts, hemangioblastoma-derived hemangioblasts demonstrated self-renewal and differentiated into hematopoietic (erythroid and myeloid) and endothelial lineages in this environment. These findings are in agreement with observations in other tumors in which stem cells have been proposed as the origin of tumor initiation and propagation [27,[38][39][40][41].
According to Knudson's ''two-hit'' hypothesis of tumorigenesis [42], initiation of tumor formation follows inactivation of both alleles of the VHL gene. Because germline mutations of the VHL gene were present in all of these patients. and deletion of the remaining wild-type allele is necessary to develop VHL-associated tumors (i.e., hemangioblastomas), we were able to confirm the origin of the hematopoietic and endothelial progeny by loss of heterozygosity analyses. Deletion of the wild-type VHL allele (loss of heterozygosity) was found in the cultured hematopoietic and endothelial progeny, confirming their neoplastic origin.
The functional ability of the hemangioblast-derived immature hematopoietic progenitors to expand into mature erythrocytic and granulocytic progeny was tested using a culture media designed to expand immature hematopoietic cells ex vivo [14] and to exploit the characterized receptor expression (erythropoietin receptor, c-kit, and IL-3). The capacity of the hemangioblast-derived endothelial progenitors to proliferate in a manner similar to early vasculogenesis was examined using a different culture medium designed to generate endothelial cells ex vivo and that also exploited the characterized receptor expression (erythropoietin receptor and Flk-1 [VEGF receptor-2]) [17]. As occurs with extramedullary hematopoiesis and blood island formation, expansion of mature enucleated erythrocytes and endothelial cells occurred in these environments. These findings are consistent with the presence of a developmentally arrested hemangioblast that can be functionally reactivated under appropriate environmental conditions and concur with histological evidence of intratumoral blood island formation in hemangioblastomas [12,13].

Implications for Clinical Development of Hemangioblastomas
While our findings indicate that the cell of origin for VHLassociated hemangioblastomas are established during development, clinically significant or radiographic-evident hemangioblastomas rarely occur in infancy or early childhood in VHL. This observation is consistent with the known tendency of these tumors for slow growth and well-documented quiescent growth periods that can last for years [11,43]. Thus, it is likely that in early childhood many of these tumors remain asymptomatic and below the level of radiographic detection (even by sensitive MRI techniques). An analogous situation is the presence of endolymphatic sac tumors, which are undetectable by imaging for years in many patients with VHL [5,6,44].
The presence of subclinical and infraradiologic hemangioblastoma tumorlets in VHL patients at autopsy supports the concept that these tumors may be developmentally arrested but reactivated under appropriate conditions later in life [34]. There are several factors that could conceivably be involved in the propagation of these tumors. Circulating factors, including hormones (e.g., during puberty or pregnancy) could promote the growth of previously transformed tumor cells [45][46][47]. Similarly, local microenvironmental stimuli driven by autocrine or paracrine loops have been hypothesized to spur hemangioblastoma growth [13]. Finally, additional genetic alterations could underlie the observed delayed appearance and erratic growth of these tumors.

Implications for the Distribution of Hemangioblastomas
Because embryologic development results in the precise topographical arrangement of various cell types, a hemangioblast-defined origin of hemangioblastomas could explain the distribution of these tumors. We recently examined the natural history and location of CNS hemangioblastomas (655 CNS hemangioblastomas in 160 consecutive patients) in VHL patients as defined by MRI [11]. Consistent with the distribution of hemangioblastomas found here, that study revealed that hemangioblastomas were found only within very defined regions of the CNS, including the cerebellum (38%), brainstem (10%), and spinal cord (51%).
This tumor distribution coincides with previous topographical analyses, which revealed that Scl is transiently expressed by embryologic hemangioblast cells during development in the retina, diencephalon, mesencepahlon, metencephalon, and spinal cord, but not the telencephalon [22,24,[48][49][50]. Thus, the highly specific distribution of Scl gene expression during embryogenesis correlates with the restricted distribution of hemangioblastomas in VHL patients [11]. The rare occurrence (less than 1%) of hemangioblastomas in the supratentorial cerebral hemispheres may represent an ectopic rest of hemangioblasts. Taken together, these findings suggest that hemangioblastoma distribution in VHL does not result from migration of the tumor-initiating cells (hemangioblast); instead, developmental processes may govern distribution. The topographic and tissue-restricted distribution patterns of tumors that occur in other familial neoplasia syndromes may also arise by a similar mechanism.

Study Limitations
Because we examined only VHL-associated hemangioblastomas, the results of this study may not be applicable to hemangioblastomas that arise sporadically. Moreover, because hemangioblastomas are benign neoplasms, these findings may not apply to malignant tumors, particularly those that do not arise in the context of tumor suppressor syndromes. Future analyses of other malignant tumors may be necessary to understand the factors involved with growth initiation, clinical development, and the developmental effects of tumor suppressor genes.

Potential Clinical Implications
Several findings of this study may have future implications for clinical practice. The data here support the concept of a multipotent tumor stem cell as the initiating and sustaining cell in hemangioblastomas in the context of VHL. Specifically, an embryologically arrested hemangioblast that retains the ability to form blood and endothelial cells appears to underlie the formation and propagation of benign VHLassociated hemangioblastomas. Identification of the hemangioblast as the cell of origin in hemangioblastomas could permit the selective targeting of uniquely expressed proteins or modulation of specific molecular pathways in the selective treatment of these neoplasms. Moreover, strategies that exploit the known differentiation potential (blood and endothelial cells) of the hemangioblast may offer an alternative treatment paradigm.

Conclusions
These findings indicate that the neoplastic cells from CNS hemangioblastomas from VHL patients are hemangioblasts derived from embryologically arrested mesoderm. These embryologic tumor cells retain their multipotent differentiation ability and can develop into hematopoietic or endothelial progeny. The identification of the cell of origin of this tumor enhances understanding of the biological basis of this familial tumor suppressor syndrome and may provide therapeutic opportunities.

Accession Numbers
The GenBank (http://www.ncbi.nlm.nih.gov) Entrez Gene identification number for VHL is 7428. coexpressed with GATA-1 in hemopoietic cells but is also expressed in developing brain.

Editors' Summary
Background. von Hippel-Lindau (VHL) disease is a rare genetic condition characterized by the development of benign and malignant tumors in multiple organ systems. All the cells of people with this disorder contain one normal copy of the VHL gene and one altered copy. This gene encodes a tumor suppressor, a protein that stops tumors growing. One functioning copy of the VHL gene is sufficient to prevent any problems, but if the remaining normal copy becomes altered (mutated) in individual cells of patients with VHL disease, tumors result. These tumors are mainly benign (noncancerous growths that do not spread around the body) and form in parts of the body that are rich in blood vessels, in particular in the retina (the back of the eye), the cerebellum (the back of the brain), the brainstem (which links the cerebellum to the spinal cord), and the spinal cord. These central nervous system (CNS) tumors are called hemangioblastomas and look like little knots of capillaries (fine blood vessels). As they grow, they can cause problems through fluid leakage or by pressing on brain tissue. There is no cure for VHL disease, but patients can be monitored and their tumors dealt with before they get too large.
Why Was This Study Done? It is not known what sort of cells hemangioblastomas develop from or why they occur only in specific parts of the CNS. This information could help researchers develop ways to prevent or treat these hemangioblastomas. One possibility is that hemangioblastomas develop from a special kind of embryonic cell called a hemangioblast. This multipotent stem cell-a continuously dividing cell that can develop (differentiate) into several nondividing cell typesis the source of blood cells and blood vessel cells in the embryo. In this study, the researchers have examined tumor cells (so-called neoplastic stromal cells) taken from hemangioblastomas to see whether this theory is correct.
What Did the Researchers Do and Find? The researchers obtained several CNS hemangioblastomas from patients with VHL disease and stained slices of them with antibodies that stick to proteins made only by specific types of embryonic cells. This experiment showed that the neoplastic stromal cells in the tumors contained two proteins (brachyury and Flk-1) that mark hemangioblast cells in embryos and a protein called Scl that is required for blood cell formation. The neoplastic stromal cells also made several proteins expressed by the precursors of different blood cell types. When the researchers grew neoplastic stromal cells from the hemangioblastomas in different conditions, they found that the cells differentiated into the precursors of two types of blood cell (erythrocytes and granulocytes) and of the cells that line blood vessels (endothelial cells). The researchers confirmed that these precursors had arisen from the neoplastic cells in the hemangioblastomas by showing that they did not contain a normal copy of the VHL gene.
What Do These Findings Mean? These findings indicate that the neoplastic stromal cells in CNS hemangioblastomas in patients with VHL disease are indeed hemangioblasts. Their expression of the protein brachyury, which is normally only expressed in early development, suggests that these hemangioblasts were arrested during embryogenesis, possibly when they lost their normal VHL gene. Because VHLrelated tumors usually appear in adulthood, circulating factors such as hormones might promote the growth of these pre-existing but dormant tumor cells at specific stages of life. Furthermore, the expression of Scl in the neoplastic stromal cells suggests that the specific pattern of hemangioblastomas in VHL disease is determined during development since Scl is transiently expressed in hemangioblasts in the parts of the developing brain that subsequently harbor hemangioblastomas. Because the researchers only examined VHL-associated hemangioblastomas, their findings do not explain the origin of hemangioblastomas in patients without VHL disease or the origin of VHL-associated malignant tumors. However, the identification of hemangioblasts as the source of hemangioblastomas in VHL disease suggests that it might be possible to treat these tumors by persuading the hemangioblasts to differentiate or by using their unique proteins to target toxic drugs to hemangioblastomas.