BACKGROUND
nearly all eukaryotic messenger rna  carries a long series of adenine at the 3' end called the polyadenylation ) tail.
the molecular process synthesizing the poly tail is called polyadenylation.
eukaryotic polyadenylation was first reported more than half a century ago  <cit> .
since then, tremendous progress has been made in elucidating the mechanism, regulation, protein factors, and related biological functions.
although polyadenylated transcripts in prokaryotes were first identified since  <dig>  <cit> , the majority of studies focus on eukaryotes and their dna viruses, probably due to the obstacles of isolating unstable prokaryotic transcripts.
more recently, polyadenylation has been studied in archaea  <cit>  and in organelles: the chloroplast  <cit> , and mitochondria  <cit> .
the prevalence of polyadenylation across all three domains of life signifies a long evolutionary history in which varied complexity and additional functions have been selected by diverse species.
polyadenylation consists of two tandem enzymatic reactions: the cleavage of a nascent mrna from the elongating rna polymerase, followed by the non-template synthesis of a poly tail that varies in length between speices.
a typical eukaryotic poly site is characterized by three cis-elements.
the first element lies where the pre-mrna is cut off from the rna polymerase at the pre-mrna's 3'-most exon: the cleavage site.
the second element is a highly conserved hexanucleotide, namely the poly signal.
the majority of poly signals are located ~ <dig> nts upstream from the cleavage sites.
66% and 16% of mammalian transcripts contain aauaaa and auuaaa, respectively  <cit> , making the canonical poly signal awuaaa .
the third element is named the downstream element  which is located at ~10- <dig> nts downstream from the cleavage site.
in contrast to the poly signal, no consensus sequence has been found in the dse among animals except that it is enriched mainly with 'u' and sprinkled with 'g'.
therefore the dse is known as u/gu-rich region.
although cis-elements are short and variable, polyadenylation takes place precisely  at the same location  of a gene.
moreover, even though all genes within a species are processed by the same set of core polyadenylation factors, two poly sites rarely resemble each other  <cit> .
the functionally conserved but sequence-variable poly sites not only challenge the identification of definitive features for recognition, but also present an intriguing case study for the understanding of the evolution of non-coding regions in different species.
we present an improved poly site model that distinguishes itself from existing models in four ways.
1) instead of choosing features haphazardly, we use principal component analysis  to identify the localization of cis-elements without presuming what they are.
2) our four feature model uses fewer features than existing methods , which use between six and over  <dig>  features  <cit> , and achieves superior prediction accuracy.
the rationale of taking a parsimonious approach in feature selection is to circumvent the dimensionality curse  <cit> , but our simple model also requires a smaller training dataset as a result.
3) despite the highly variable poly site cis-elements, the poly complex is still able to cleave the transcript at the same position.
we believe the poly site is marked by more information than just sequence elements, such as peculiar chromatin structure  <cit> .
therefore, we have incorporated nucleosome occupancy as a novel feature in our model.
4) we have used seven diverse species to validate the generality of our four-feature model, a far wider range of species than has been attempted when validating existing methods.
the seven species are homo sapiens , mus musculus , gallus gallus , caenorhabditis elegans , oryza sativa , arabidopsis thaliana , and solanum lycopersicum .
intriguingly, the performance of our model on cross species predictions reflected the phylogenetic distances among these seven eukaryotes.
materials and methods
notation of region
we use the notation < ± n, ± m > to denote a region with respect to the cleavage site, where n and m are the number of nucleotides  upstream  or downstream , from the cleavage site.
poly site discovery
poly sites were discovered by mapping polyadenylated ests and/or cdnas to the reference genomes as described in  <cit> .
briefly, our method involves: a) ests that either terminate with at least eight polymeric 'a' at the 3' end, or start with at least eight polymeric 't' at the 5' end are selected.
for cdna datasets, only sequences that terminate with at least eight 'a' are selected.
b) polyadenylated ests and cdnas are mapped to the appropriate reference genome using ncbi blast  <dig> .
<dig>  c) customized python scripts are used to determine the direction of transcription, and to eliminate artifacts due to false oligo-dt priming.
ests and cdnas sequences were downloaded from ncbi's dbest, and refseq databases.
through this method,  <dig> ,  <dig> ,  <dig> ,  <dig>  and  <dig>  poly sequences were discovered in human, mouse, chicken, c.elegans, and tomato, respectively.
genomes of human build  <dig>  mouse build  <dig>  chicken may  <dig> release, c.elegans feb  <dig> release, and tomato nov  <dig> release  <dig>  were used for mapping.
other poly sequences
we used reliable published datasets of arabidopsis poly sequences from http://www.users.muohio.edu/liq/loke_et_al_2005_8k <cit>  and rice from http://www.users.muohio.edu/liq/rice_55k_polya_site_dataset.zip <cit> .
we eliminated poly sequences that could not be mapped to the reference genomes.
arabidopsis genome build  <dig>  <dig> and rice genome build  <dig>  jun  <dig> were used.
as a result, we compiled  <dig>  and  <dig>  poly sequences from arabidopsis and rice, respectively.
transcribed non-poly sequences
for specificity testing,  <dig>  spliced and unspliced transcribed sequences were chosen randomly from arabidopsis and human datasets.
gene sequences of arabidopsis were downloaded from ftp://ftp.arabidopsis.org/home/tair/sequences/blast_datasets/tair10_blastsets/, whereas human sequences were prepared by customized python programs.
for sequences longer than  <dig> nts, 600-nt long fragments were extracted at random locations.
position-by-kmer matrix
in order to examine the localization of kmers  in a set of sequences, we captured and converted sequence information into a matrix so that it could be processed by pca .
starting from the leftmost position of a sequence, we use a sliding window of size k to count the occurrence of a kmer and its position.
the occurrence of each kmer at various positions is stored in the column of the matrix, so that each row stores the counts of all possible kmers  at a position.
a row is a 4k -dimensional vector, where each dimension is the occurrence of a kmer at a position.
each column is the occurrences of a kmer at different positions.
position score matrix
position score matrix  is a m-by-n matrix, where m is equal to 4k, k is the size of oligonucleotide.
and, n is , where l is the length of considered region.
for our model construction, we consider region < - <dig>  + <dig> >, i.e.
l is equal to  <dig>  the value of each entry in psm, v, is defined as below:
v=log2observedrealkmer_at_position_nobservedfalsekmer_at_position_n
the value v represents the fold difference of observing a particular kmer at position n in real versus false poly sequences, where n is in the range of  <dig> to .
given a sequence, a k-sized window is used to slide through it one position at a time from left to right.
the result is a list of overlapping kmers with their corresponding starting positions in the sequence  pairs).
this list will be used to retrieve the associated values from psm i.e.
v, and the sum of these values is the psm score of the sequence.
nucleosome occupancy matrix
nom is constructed by a similar method used in building the psm.
the nucleosome prediction method developed by segal lab  <cit>  was used to determine the probability of nucleosome occupancy in each position in < - <dig>  + <dig> > for real and false poly sequences.
due to boundary effects, predictions from  <dig> nts at both ends are discarded, leaving only the middle < - <dig>  + <dig> > region for consideration.
like psm, the value of each entry in the nom, v, is a defined as below
v=log2occupancy_prob_at_position_nrealoccupancy_prob_at_position_nfalse
pn is the predicted occupancy probability at position n.
feature vector
the feature vector consists of four features: trimer psm scores for < - <dig>  - <dig> > and < + <dig>  + <dig> >, dimer psm score across the cleavage site < - <dig>  + <dig> >, and the nucleosome occupancy score based on nom in the region < - <dig>  + <dig> >.
feature vectors are standardized by the mean and standard deviation obtained from all real and false poly sequences within each species' dataset.
performance measurement
prediction is measured by sensitivity , specificity , and matthews correlation coefficient  <cit> .
sn = tp/, sp = tn/, mcc = /√•••, where tp, tn, fp, and fn denote true positive, true negative, false positive, and false negative, respectively.
training and testing procedures
false poly sequences  were generated according to the 2nd order markov model obtained from real poly sequences in region < - <dig>  + <dig> >.
accuracy was reported as the average of twenty ten-fold cross validations.
linear discriminant analysis  and logistic regression  were performed by python machine learning package mlpy  <cit> .
method comparisons
polya_svm  <cit>  version  <dig>  was downloaded from http://exon.umdnj.edu/polya_svm/, and polya-ep  <cit>  was downloaded from http://mlkd.csd.auth.gr/polya/tools.html.
polya_svm was tested on  <dig>  human poly sequences, each consisted of  <dig> nts upstream and downstream from the poly site.
for polya-ep,  <dig>  arabidopsis poly sequences were used where each consisted of the  <dig> nts upstream and  <dig> nts downstream regions, relative to cleavage sites.
18s, gapdh, and cpsf <dig> sequences, and phylogenetic distance calculations
the evolutionary relationship among these seven eukaryotes was determined through phylogenetic analysis of three independent gene products: the ribosomal 18s, gylceraldehyde  <dig> phosphate dehydrogenase  and cleavage and polyadenylation specificity factor subunit  <dig> .
18s genes of human , mouse , c. elegans , o. sativa , s. lycopersicum , and a. thaliana  were obtained from ncbi.
chicken's 18s was determined by a blat search  <cit>  using human's 18s as query sequence in ucsc genome browser, http://genome.ucsc.edu <cit> .
homologous gapdh protein sequences for all organisms except s. lycopersicum were obtained from the homologene database in ncbi .
the gapdh sequence of s. lycopersicum was obtained by translating its cdna sequence ak <dig>  similarly, homologous cpsf <dig> protein sequences  were collected from the homologene database, except for s.lycopersicum, where the cpsf <dig> was obtained by the translation of ak <dig>  sequences were aligned by t_coffee  <cit> .
phylogenetic distances between the seven species were calculated by dnadist  and protdist  in phylip  <cit> .
RESULTS
feature identification by principal component analysis
simulated sequences were used to prove the viability of our pca method , which we then applied to arabidopsis and human data using different sizes of kmers.
figure 1a shows pca of arabidopsis poly sequences using dimers.
in the pca-position profile , four localized groups  are identified: - <dig> , - <dig> , - <dig> , and + <dig> .
examining the pca-oligo profile , 'aa' and 'at' are most often found in the dur.
'tt' is the most significant dimer in the pur and dr. note that the extent of localization in dr is moderate versus far upstream and downstream regions.
this result concurs with current view that poly sites of plants are mainly characterized by upstream elements  <cit> .
'ca', 'ga', and 'ta'  are found in cs, and of these, 'ca' is predominant.
similar to arabidopsis, human  pur is also found to localize in four regions, but with different relative concentration of dimers.
dur at - <dig> shows the strongest localization of 'aa', and then followed by 'at' and 'ta'.
pur at - <dig> shows a weak localization  of 'tt'.
the human cs at + <dig> shows weaker localization compared with arabidopsis' cs.
however, seven dimers  are likely to localize at the human cs .
the human dr covers an extensive region from + <dig> to + <dig>  which is absent in arabidopsis.
'ct', 'gt', 'tc', 'tg', and 'tt' are the main dimers in human dr.
even with trimer motifs, there are no prominent cis-elements in the downstream region of arabidopsis poly sites.
the downstream elements  are not essential for polyadenylation in the invertebrate c.elegans  <cit> .
these results might indicate that dse, which are essential in vertebrates, evolved more recently than other poly site elements.
moreover, we do not find core poly cis-elements stretching beyond  <dig> nts upstream from the cleavage site in arabidopsis, which is at odds with an early bioinformatics study  <cit> .
certain genes do possess specific cis-elements in a much further upstream region  <cit> , but there weren't sufficient signals across all genes for this to be a key feature of poly sites.
in order to construct a model that can embrace such a diversity in the poly sites in different species, the features obtained from the region flanking the poly site are split into three subregions approximately according to the peaks as shown in figure  <dig>  i.e.
< - <dig>  - <dig> >, < - <dig>  + <dig> >, and < + <dig>  + <dig> >.
the purpose of doing that is to prevent the contributing effect of features from one subregion to be cancelled by non-contributing features from other subregions.
nevertheless, poly sites may be constituted of longer sequence structures, and we continued pca for longer kmers.
differences between dimer and trimer profiles  are insignificant in both arabidopsis and humans.
however, when the size of kmers is increased to six, the localization signal at the cs drops or even vanishes in both species , meaning that the size of elements favored at the cs is shorter than  <dig> nts, and may be identified well by dimers or trimers.
with hexamers, the significant signal in the human pur vanishes  but not for arabidopsis.
the dr of human becomes flat relative to the dur and its signal disappears entirely in the octamer pca-position profile.
hence, the size of sequence elements in the human dr should be 3- <dig> nts.
in both species, the canonical poly signal aauaaa has the longest pca-distance according to the pca-oligo profiles in figure s2b.
interestingly, unlike humans, the dur and pur of arabidopsis still remain even in octamer pca-position profile .
as the poly signal is highly conserved, the hexamer profiles presented above may lead one to believe that the poly signal and/or hexamers are effective features for poly site recognition.
on the contrary, our study suggests the opposite .
the matthews correlation coefficient  attained by the trimer model was  <dig>  versus  <dig>  for the hexamer model, indicating the trimer model had more accurate predictive power.
furthermore, we compared the false positive rates of predicting false poly sequences between the trimer model and an existing method, polya_svm  <cit> , which uses a hexamer as an upstream feature.
surprisingly, the false positive rate for the hexamer feature-containing polya_svm was 82%, while that of the trimer model was just 17%.
clearly, hexamer models do not offer substantial advantages relative to a trimer model.
additionally, the model with hexamer features demands a larger training sample  than the trimer model because smaller samples are unlikely to capture the full distribution of all genomic contexts for poly hexamers within a species.
complicating hexamer model training is the fact that poly signallike hexamers are ubiquitous in large genomes in places that are not poly sites, such as introns.
eukaryotic cells likely enhance the recognition and regulation of the actual poly sites with other auxiliary cis-elements that may be shorter than hexamers  <cit> , which could be better distinguished with a model looking for smaller features.
this viewpoint is further supported by the fact that intra-species poly sites are highly variable and yet they are recognized by the same polyadenylation apparatus.
thus, we choose trimers instead of hexamers to form part of the feature vector.
nucleosome structure at poly sites
previous studies suggested that nucleosomes are depleted in the proximity of poly sites  <cit> , but no poly site prediction models have yet incorporated this feature.
hence we utilize a nucleosome occupancy prediction method developed by the segal lab  <cit>  to explore nucleosome occupancy in the region flanking the poly sites .
even though the false poly sequences mimic the 2nd order markov property of real poly sequences, their predicted nucleosome occupancy remains at a steady level throughout the entire region.
in contrast, both human and arabidopsis actual poly sequences show a reduced likelihood of nucleosomes around the cleavage site.
our results indicate a putative relationship between nucleosome depletion and polyadenylation site, in agreement with published works  <cit> .
we have thus included predicted occupancy in the feature vector.
the role of nucleosome formation on polyadenylation is currently unknown, and should be further investigated.
we speculate nucleosome formation may influence the selection of alternative poly sites in genes, as 54%, 32% and 70% of human, mouse, and arabidopsis genes, respectively, contain multiple poly sites  <cit> .
feature vector and poly site predictions
after a thorough analysis of poly sequences, we have selected four features in the feature vector to signify a poly site .
we validated the efficacy of these features by two independent supervised machine learning  methods: logistic regression  and linear discriminant analysis .
lda requires a gaussian distribution of feature values but lr does not.
quadratic discriminant analysis  was also considered since the positive and negative datasets might exhibit different covariance, but we did not find any significant differences between lda and qda in terms of prediction accuracy .
for simplicity, we report our results using lda.
the optimal performance peaks when thresholds are between  <dig>  and  <dig>  according to the receiver operating characteristic  method , therefore the threshold for lr is set to  <dig>  for the rest of the analysis.
the two ml methods were tested using poly sequences from seven species including two mammals, a bird, an invertebrate, two dicots, and a monocot.
our two ml methods achieved 85-92% specificity and 85-96% sensitivity for each of the seven species .
results from the logistic and lda show no significant differences .
since no existing method has been used for multispecies predictions, we must compare our methods to other methods on individual species.
sn and sp are expressed in percentage.
methods comparison and model validation
we compared our lr and lda models to two poly site prediction methods, polya_svm  <cit> , and polya-ep  <cit> .
the former was originally developed for animals and the latter was for plants.
both of our simpler, four-feature lda and lr models produced better results than polya_svm in predicting human poly sites in terms of sensitivity and specificity .
the polya-ep model is skewed towards sensitivity when analyzing arabidopsis data, while our models performed well and without bias.
as our models were trained on the arabidopsis dataset built by the pac method  <cit> , we can also compare our two ml methods with this third published method without repeating the tests.
the best mcc achieved by pac is between  <dig>  to  <dig>  with various types of negative samples and positions of deviation .
but our two ml methods are able to improve mcc significantly to  <dig> - <dig>  .
additionally, we conducted three more tests to validate our model: i) prediction for other transcribed, non-poly genomic sequences, ii) testing the relative contribution of each feature to accuracy, and iii) prediction of larger poly regions.
these results can be found in the additional file  <dig>
sn and sp are expressed in percentages.
the mcc of pac is from figure  <dig> of  <cit> .
cross species poly site predictions
as polyadenylation is a universal function across the tree of life, we are interested in the evolution of poly sites across diverse species.
however, because poly sites lie in relatively unconstrained untranslated regions , aligning these sequences is usually infeasible.
as a result, prevalent phylogenetic methods are unsuitable to look at the evolution of poly sites.
our method, applied across species, is able to recapitulate the phylogenetic relationship across species based on non-coding sequences.
the phylogenetic distances among these seven organisms  was established through phylogenetic trees, built independently using three unlinked genes .
all three phylogenetic trees not only produced congruent topology, but also similar evolutionary distances .
unit is substitution per site.
the last column rsn is defined according to materials and methods section.
we utilize the reciprocal sensitivity   to assess the overall homology of poly sites between species.
sensitivity was obtained by applying the logistic model of one species to predict real poly site sequences of another species .
the validity of the poly-based phylogenetic relationship is measured by the correlation between gene-based phylogenetic distances and rsn predicted by the logistic model.
the correlation coefficients are - <dig> , - <dig> , and - <dig>  based on the 18s, gapdh, and cpsf <dig> trees, respectively  .
given the high intraspecies variability, nearly no two poly sites are alike in the same genome  <cit> , there should therefore not necessarily be conservation in poly sites between more closely related species.
however, our results show that prediction accuracy decreases when the model from one species is applied to data from other species and the decrease is proportional to phylogenetic distance.
sensitivity is calculated by using model from one species  to predict real and false poly sequences from other species .
the table needs not be symmetrical because poly sites from different species tend to possess different characteristics according to the pca profiles discussed above.
this leads us to question about the boundary of potential poly sites, which we can conceptualize as pa-space, the allowable variations of functional poly sites.
core polyadenylation proteins that directly bind to the pre-mrna, such as cpsf <dig> and cstf- <dig>  are likely to co-evolve with the pa-space.
while some studies have identified conserved residues in cstf- <dig> that correlated with poly downstream elements  <cit> , no one has looked at how variation in the nucleotide sequence and proteins are correlated and coevolving.
currently, the sequence-level dynamics between rna binding proteins and their substrates is unclear, but the decreasing costs of sequencing transcriptomes should provide data from a wider range of species soon.
CONCLUSIONS
we have shown the feasibility of conducing comprehensive genomic analysis of poly sites using pca, a method which could be broadly applied for any cis-element identification.
we believe a model focused on very short oligonucleotides outperformed those with hexamers as features because it embraces the conspicuous poly signal elements without sacrificing the diverse family of auxiliary biological signals surrounding the poly sites.
we also included for the first time nucleosome occupancy as an informative predictor of poly sites.
list of abbreviations
3'utr: 3' untranslated region; cs: cleavage site; dr: downstream region; dse: downstream element; dur: distal upstream region; fn: false negative; fp: false positiv;.
lda: linear discriminant analysis; lr: logistic regression; mcc: matthews correlation coefficien;.
ml: machine learning; nom: nucleosome occupancy matrix; pca: principal component analysis; poly: polyadenylation; psm: position score matrix; pur: proximal upstream regio; qda: quadratic discriminant analysis; roc: receiver operating curve; rsn: reciprocal sensitivity; sn: sensitivity; sp: specificity; tn: true negative; tp: true positive.
competing interests
the authors declare that they have no competing interests.
authors' contributions
esh conceived of the study, designed and carried out the analyses and drafted the manuscript.
sig conceived of the study.
sd conceived of the plant portion of the study, assisted in data analysis and helped draft the manuscript.
supplementary material
additional file 1
table s1: existing poly sites methods.
figure s1: feature identification by principal component analysis .
figure s2: pca of real poly sequences.
the canonical poly signal is a misguiding feature.
figure s3: roc of logistic method.
table s2: false positive rate  committed by different methods in handling cds gene sequences.
table s3: relative contribution of individual features.
figure s4: predictions of five hundred  <dig> -nt poly sites by sliding a 600-nt sliding window from left to right.
figure s <dig>  correlation between phylogenetic distance and reciprocal sensitivity  between seven species: human, mouse, chicken, c.elegans, rice, arabidopsis, and tomato.
click here for file
declarations
the publication fee of this article is funded by rutgers school of environmental and biological sciences and the new jersey agricultural experiment station .
this work is supported by national institute of health  grant nih k <dig> gm093854- <dig> to e.s.h.
s.i.g.
is supported by nih gm <dig>  s.d.
is supported by funds from the njaes.
this article has been published as part of bmc bioinformatics volume  <dig> supplement  <dig>  2013: selected articles from the eleventh asia pacific bioinformatics conference : bioinformatics.
the full contents of the supplement are available online at http://www.biomedcentral.com/bmcbioinformatics/supplements/14/s <dig>

 nearly all eukaryotic messenger rna carries a long series of adenine at the 3' end called the polyadenylation ) tail. the molecular process synthesizing the poly tail is called polyadenylation. eukaryotic polyadenylation was first reported more than half a century ago. since then, tremendous progress has been made in elucidating the mechanism, regulation, protein factors, and related biological functions. although polyadenylated transcripts in prokaryotes were first identified since, the majority of studies focus on eukaryotes and their dna viruses, probably due to the obstacles of isolating unstable prokaryotic transcripts. more recently, polyadenylation has been studied in archaea and in organelles: the chloroplast, and mitochondria. the prevalence of polyadenylation across all three domains of life signifies a long evolutionary history in which varied complexity and additional functions have been selected by diverse species. polyadenylation consists of two tandem enzymatic reactions: the cleavage of a nascent mrna from the elongating rna polymerase, followed by the non-template synthesis of a poly tail that varies in length between speices. a typical eukaryotic poly site is characterized by three cis-elements. the first element lies where the pre-mrna is cut off from the rna polymerase at the pre-mrna's 3'-most exon: the cleavage site. the second element is a highly conserved hexanucleotide, namely the poly signal. the majority of poly signals are located ~ nts upstream from the cleavage sites. 66% and 16% of mammalian transcripts contain aauaaa and auuaaa, respectively, making the canonical poly signal awuaaa. the third element is named the downstream element which is located at ~10- nts downstream from the cleavage site. in contrast to the poly signal, no consensus sequence has been found in the dse among animals except that it is enriched mainly with'u' and sprinkled with'g'. therefore the dse is known as u/gu-rich region. although cis-elements are short and variable, polyadenylation takes place precisely at the same location of a gene. moreover, even though all genes within a species are processed by the same set of core polyadenylation factors, two poly sites rarely resemble each other. the functionally conserved but sequence-variable poly sites not only challenge the identification of definitive features for recognition, but also present an intriguing case study for the understanding of the evolution of non-coding regions in different species. we present an improved poly site model that distinguishes itself from existing models in four ways. 1) instead of choosing features haphazardly, we use principal component analysis to identify the localization of cis-elements without presuming what they are. 2) our four feature model uses fewer features than existing methods, which use between six and over features, and achieves superior prediction accuracy. the rationale of taking a parsimonious approach in feature selection is to circumvent the dimensionality curse, but our simple model also requires a smaller training dataset as a result. 3) despite the highly variable poly site cis-elements, the poly complex is still able to cleave the transcript at the same position. we believe the poly site is marked by more information than just sequence elements, such as peculiar chromatin structure. therefore, we have incorporated nucleosome occupancy as a novel feature in our model. 4) we have used seven diverse species to validate the generality of our four-feature model, a far wider range of species than has been attempted when validating existing methods. the seven species are homo sapiens, mus musculus, gallus gallus, caenorhabditis elegans, oryza sativa, arabidopsis thaliana, and solanum lycopersicum. intriguingly, the performance of our model on cross species predictions reflected the phylogenetic distances among these seven eukaryotes. materials and methods notation of region we use the notation < ± n, ± m > to denote a region with respect to the cleavage site, where n and m are the number of nucleotides upstream or downstream, from the cleavage site. poly site discovery poly sites were discovered by mapping polyadenylated ests and/or cdnas to the reference genomes as described in. briefly, our method involves: a) ests that either terminate with at least eight polymeric'a' at the 3' end, or start with at least eight polymeric't' at the 5' end are selected. for cdna datasets, only sequences that terminate with at least eight'a' are selected. b) polyadenylated ests and cdnas are mapped to the appropriate reference genome using ncbi blast. c) customized python scripts are used to determine the direction of transcription, and to eliminate artifacts due to false oligo-dt priming. ests and cdnas sequences were downloaded from ncbi's dbest, and refseq databases. through this method,,,, and poly sequences were discovered in human, mouse, chicken, c.elegans, and tomato, respectively. genomes of human build mouse build chicken may release, c.elegans feb release, and tomato nov release were used for mapping. other poly sequences we used reliable published datasets of arabidopsis poly sequences from http://www.users.muohio.edu/liq/loke_et_al_2005_8k and rice from http://www.users.muohio.edu/liq/rice_55k_polya_site_dataset.zip. we eliminated poly sequences that could not be mapped to the reference genomes. arabidopsis genome build and rice genome build jun were used. as a result, we compiled and poly sequences from arabidopsis and rice, respectively. transcribed non-poly sequences for specificity testing, spliced and unspliced transcribed sequences were chosen randomly from arabidopsis and human datasets. gene sequences of arabidopsis were downloaded from ftp://ftp.arabidopsis.org/home/tair/sequences/blast_datasets/tair10_blastsets/, whereas human sequences were prepared by customized python programs. for sequences longer than nts, 600-nt long fragments were extracted at random locations. position-by-kmer matrix in order to examine the localization of kmers in a set of sequences, we captured and converted sequence information into a matrix so that it could be processed by pca. starting from the leftmost position of a sequence, we use a sliding window of size k to count the occurrence of a kmer and its position. the occurrence of each kmer at various positions is stored in the column of the matrix, so that each row stores the counts of all possible kmers at a position. a row is a 4k -dimensional vector, where each dimension is the occurrence of a kmer at a position. each column is the occurrences of a kmer at different positions. position score matrix position score matrix is a m-by-n matrix, where m is equal to 4k, k is the size of oligonucleotide. and, n is, where l is the length of considered region. for our model construction, we consider region < - + >, i.e. l is equal to the value of each entry in psm, v, is defined as below: v=log2observedrealkmer_at_position_nobservedfalsekmer_at_position_n the value v represents the fold difference of observing a particular kmer at position n in real versus false poly sequences, where n is in the range of to. given a sequence, a k-sized window is used to slide through it one position at a time from left to right. the result is a list of overlapping kmers with their corresponding starting positions in the sequence pairs). this list will be used to retrieve the associated values from psm i.e. v, and the sum of these values is the psm score of the sequence. nucleosome occupancy matrix nom is constructed by a similar method used in building the psm. the nucleosome prediction method developed by segal lab was used to determine the probability of nucleosome occupancy in each position in < - + > for real and false poly sequences. due to boundary effects, predictions from nts at both ends are discarded, leaving only the middle < - + > region for consideration. like psm, the value of each entry in the nom, v, is a defined as below v=log2occupancy_prob_at_position_nrealoccupancy_prob_at_position_nfalse pn is the predicted occupancy probability at position n. feature vector the feature vector consists of four features: trimer psm scores for < - - > and < + + >, dimer psm score across the cleavage site < - + >, and the nucleosome occupancy score based on nom in the region < - + >. feature vectors are standardized by the mean and standard deviation obtained from all real and false poly sequences within each species' dataset. performance measurement prediction is measured by sensitivity, specificity, and matthews correlation coefficient. sn = tp/, sp = tn/, mcc = /√•••, where tp, tn, fp, and fn denote true positive, true negative, false positive, and false negative, respectively. training and testing procedures false poly sequences were generated according to the 2nd order markov model obtained from real poly sequences in region < - + >. accuracy was reported as the average of twenty ten-fold cross validations. linear discriminant analysis and logistic regression were performed by python machine learning package mlpy. method comparisons polya_svm version was downloaded from http://exon.umdnj.edu/polya_svm/, and polya-ep was downloaded from http://mlkd.csd.auth.gr/polya/tools.html. polya_svm was tested on human poly sequences, each consisted of nts upstream and downstream from the poly site. for polya-ep, arabidopsis poly sequences were used where each consisted of the nts upstream and nts downstream regions, relative to cleavage sites. 18s, gapdh, and cpsf sequences, and phylogenetic distance calculations the evolutionary relationship among these seven eukaryotes was determined through phylogenetic analysis of three independent gene products: the ribosomal 18s, gylceraldehyde phosphate dehydrogenase and cleavage and polyadenylation specificity factor subunit. 18s genes of human, mouse, c. elegans, o. sativa, s. lycopersicum, and a. thaliana were obtained from ncbi. chicken's 18s was determined by a blat search using human's 18s as query sequence in ucsc genome browser, http://genome.ucsc.edu. homologous gapdh protein sequences for all organisms except s. lycopersicum were obtained from the homologene database in ncbi. the gapdh sequence of s. lycopersicum was obtained by translating its cdna sequence ak similarly, homologous cpsf protein sequences were collected from the homologene database, except for s.lycopersicum, where the cpsf was obtained by the translation of ak sequences were aligned by t_coffee. phylogenetic distances between the seven species were calculated by dnadist and protdist in phylip.  feature identification by principal component analysis simulated sequences were used to prove the viability of our pca method, which we then applied to arabidopsis and human data using different sizes of kmers. figure 1a shows pca of arabidopsis poly sequences using dimers. in the pca-position profile, four localized groups are identified: -, -, -, and +. examining the pca-oligo profile,'aa' and'at' are most often found in the dur.'tt' is the most significant dimer in the pur and dr. note that the extent of localization in dr is moderate versus far upstream and downstream regions. this result concurs with current view that poly sites of plants are mainly characterized by upstream elements.'ca','ga', and'ta' are found in cs, and of these,'ca' is predominant. similar to arabidopsis, human pur is also found to localize in four regions, but with different relative concentration of dimers. dur at - shows the strongest localization of'aa', and then followed by'at' and'ta'. pur at - shows a weak localization of'tt'. the human cs at + shows weaker localization compared with arabidopsis' cs. however, seven dimers are likely to localize at the human cs. the human dr covers an extensive region from + to + which is absent in arabidopsis.'ct','gt','tc','tg', and'tt' are the main dimers in human dr. even with trimer motifs, there are no prominent cis-elements in the downstream region of arabidopsis poly sites. the downstream elements are not essential for polyadenylation in the invertebrate c.elegans. these results might indicate that dse, which are essential in vertebrates, evolved more recently than other poly site elements. moreover, we do not find core poly cis-elements stretching beyond nts upstream from the cleavage site in arabidopsis, which is at odds with an early bioinformatics study. certain genes do possess specific cis-elements in a much further upstream region, but there weren't sufficient signals across all genes for this to be a key feature of poly sites. in order to construct a model that can embrace such a diversity in the poly sites in different species, the features obtained from the region flanking the poly site are split into three subregions approximately according to the peaks as shown in figure i.e. < - - >, < - + >, and < + + >. the purpose of doing that is to prevent the contributing effect of features from one subregion to be cancelled by non-contributing features from other subregions. nevertheless, poly sites may be constituted of longer sequence structures, and we continued pca for longer kmers. differences between dimer and trimer profiles are insignificant in both arabidopsis and humans. however, when the size of kmers is increased to six, the localization signal at the cs drops or even vanishes in both species, meaning that the size of elements favored at the cs is shorter than nts, and may be identified well by dimers or trimers. with hexamers, the significant signal in the human pur vanishes but not for arabidopsis. the dr of human becomes flat relative to the dur and its signal disappears entirely in the octamer pca-position profile. hence, the size of sequence elements in the human dr should be 3- nts. in both species, the canonical poly signal aauaaa has the longest pca-distance according to the pca-oligo profiles in figure s2b. interestingly, unlike humans, the dur and pur of arabidopsis still remain even in octamer pca-position profile. as the poly signal is highly conserved, the hexamer profiles presented above may lead one to believe that the poly signal and/or hexamers are effective features for poly site recognition. on the contrary, our study suggests the opposite. the matthews correlation coefficient attained by the trimer model was versus for the hexamer model, indicating the trimer model had more accurate predictive power. furthermore, we compared the false positive rates of predicting false poly sequences between the trimer model and an existing method, polya_svm, which uses a hexamer as an upstream feature. surprisingly, the false positive rate for the hexamer feature-containing polya_svm was 82%, while that of the trimer model was just 17%. clearly, hexamer models do not offer substantial advantages relative to a trimer model. additionally, the model with hexamer features demands a larger training sample than the trimer model because smaller samples are unlikely to capture the full distribution of all genomic contexts for poly hexamers within a species. complicating hexamer model training is the fact that poly signallike hexamers are ubiquitous in large genomes in places that are not poly sites, such as introns. eukaryotic cells likely enhance the recognition and regulation of the actual poly sites with other auxiliary cis-elements that may be shorter than hexamers, which could be better distinguished with a model looking for smaller features. this viewpoint is further supported by the fact that intra-species poly sites are highly variable and yet they are recognized by the same polyadenylation apparatus. thus, we choose trimers instead of hexamers to form part of the feature vector. nucleosome structure at poly sites previous studies suggested that nucleosomes are depleted in the proximity of poly sites, but no poly site prediction models have yet incorporated this feature. hence we utilize a nucleosome occupancy prediction method developed by the segal lab to explore nucleosome occupancy in the region flanking the poly sites. even though the false poly sequences mimic the 2nd order markov property of real poly sequences, their predicted nucleosome occupancy remains at a steady level throughout the entire region. in contrast, both human and arabidopsis actual poly sequences show a reduced likelihood of nucleosomes around the cleavage site. our results indicate a putative relationship between nucleosome depletion and polyadenylation site, in agreement with published works. we have thus included predicted occupancy in the feature vector. the role of nucleosome formation on polyadenylation is currently unknown, and should be further investigated. we speculate nucleosome formation may influence the selection of alternative poly sites in genes, as 54%, 32% and 70% of human, mouse, and arabidopsis genes, respectively, contain multiple poly sites. feature vector and poly site predictions after a thorough analysis of poly sequences, we have selected four features in the feature vector to signify a poly site. we validated the efficacy of these features by two independent supervised machine learning methods: logistic regression and linear discriminant analysis. lda requires a gaussian distribution of feature values but lr does not. quadratic discriminant analysis was also considered since the positive and negative datasets might exhibit different covariance, but we did not find any significant differences between lda and qda in terms of prediction accuracy. for simplicity, we report our results using lda. the optimal performance peaks when thresholds are between and according to the receiver operating characteristic method, therefore the threshold for lr is set to for the rest of the analysis. the two ml methods were tested using poly sequences from seven species including two mammals, a bird, an invertebrate, two dicots, and a monocot. our two ml methods achieved 85-92% specificity and 85-96% sensitivity for each of the seven species. results from the logistic and lda show no significant differences. since no existing method has been used for multispecies predictions, we must compare our methods to other methods on individual species. sn and sp are expressed in percentage. methods comparison and model validation we compared our lr and lda models to two poly site prediction methods, polya_svm, and polya-ep. the former was originally developed for animals and the latter was for plants. both of our simpler, four-feature lda and lr models produced better results than polya_svm in predicting human poly sites in terms of sensitivity and specificity. the polya-ep model is skewed towards sensitivity when analyzing arabidopsis data, while our models performed well and without bias. as our models were trained on the arabidopsis dataset built by the pac method, we can also compare our two ml methods with this third published method without repeating the tests. the best mcc achieved by pac is between to with various types of negative samples and positions of deviation. but our two ml methods are able to improve mcc significantly to -. additionally, we conducted three more tests to validate our model: i) prediction for other transcribed, non-poly genomic sequences, ii) testing the relative contribution of each feature to accuracy, and iii) prediction of larger poly regions. these results can be found in the additional file sn and sp are expressed in percentages. the mcc of pac is from figure of. cross species poly site predictions as polyadenylation is a universal function across the tree of life, we are interested in the evolution of poly sites across diverse species. however, because poly sites lie in relatively unconstrained untranslated regions, aligning these sequences is usually infeasible. as a result, prevalent phylogenetic methods are unsuitable to look at the evolution of poly sites. our method, applied across species, is able to recapitulate the phylogenetic relationship across species based on non-coding sequences. the phylogenetic distances among these seven organisms was established through phylogenetic trees, built independently using three unlinked genes. all three phylogenetic trees not only produced congruent topology, but also similar evolutionary distances. unit is substitution per site. the last column rsn is defined according to materials and methods section. we utilize the reciprocal sensitivity to assess the overall homology of poly sites between species. sensitivity was obtained by applying the logistic model of one species to predict real poly site sequences of another species. the validity of the poly-based phylogenetic relationship is measured by the correlation between gene-based phylogenetic distances and rsn predicted by the logistic model. the correlation coefficients are -, -, and - based on the 18s, gapdh, and cpsf trees, respectively. given the high intraspecies variability, nearly no two poly sites are alike in the same genome, there should therefore not necessarily be conservation in poly sites between more closely related species. however, our results show that prediction accuracy decreases when the model from one species is applied to data from other species and the decrease is proportional to phylogenetic distance. sensitivity is calculated by using model from one species to predict real and false poly sequences from other species. the table needs not be symmetrical because poly sites from different species tend to possess different characteristics according to the pca profiles discussed above. this leads us to question about the boundary of potential poly sites, which we can conceptualize as pa-space, the allowable variations of functional poly sites. core polyadenylation proteins that directly bind to the pre-mrna, such as cpsf and cstf- are likely to co-evolve with the pa-space. while some studies have identified conserved residues in cstf- that correlated with poly downstream elements, no one has looked at how variation in the nucleotide sequence and proteins are correlated and coevolving. currently, the sequence-level dynamics between rna binding proteins and their substrates is unclear, but the decreasing costs of sequencing transcriptomes should provide data from a wider range of species soon.  we have shown the feasibility of conducing comprehensive genomic analysis of poly sites using pca, a method which could be broadly applied for any cis-element identification. we believe a model focused on very short oligonucleotides outperformed those with hexamers as features because it embraces the conspicuous poly signal elements without sacrificing the diverse family of auxiliary biological signals surrounding the poly sites. we also included for the first time nucleosome occupancy as an informative predictor of poly sites. list of abbreviations 3'utr: 3' untranslated region; cs: cleavage site; dr: downstream region; dse: downstream element; dur: distal upstream region; fn: false negative; fp: false positiv;. lda: linear discriminant analysis; lr: logistic regression; mcc: matthews correlation coefficien;. ml: machine learning; nom: nucleosome occupancy matrix; pca: principal component analysis; poly: polyadenylation; psm: position score matrix; pur: proximal upstream regio; qda: quadratic discriminant analysis; roc: receiver operating curve; rsn: reciprocal sensitivity; sn: sensitivity; sp: specificity; tn: true negative; tp: true positive. competing interests the authors declare that they have no competing interests. authors' contributions esh conceived of the study, designed and carried out the analyses and drafted the manuscript. sig conceived of the study. sd conceived of the plant portion of the study, assisted in data analysis and helped draft the manuscript. supplementary material additional file 1 table s1: existing poly sites methods. figure s1: feature identification by principal component analysis. figure s2: pca of real poly sequences. the canonical poly signal is a misguiding feature. figure s3: roc of logistic method. table s2: false positive rate committed by different methods in handling cds gene sequences. table s3: relative contribution of individual features. figure s4: predictions of five hundred -nt poly sites by sliding a 600-nt sliding window from left to right. figure s correlation between phylogenetic distance and reciprocal sensitivity between seven species: human, mouse, chicken, c.elegans, rice, arabidopsis, and tomato. click here for file declarations the publication fee of this article is funded by rutgers school of environmental and biological sciences and the new jersey agricultural experiment station. this work is supported by national institute of health grant nih k gm093854- to e.s.h. s.i.g. is supported by nih gm s.d. is supported by funds from the njaes. this article has been published as part of bmc bioinformatics volume supplement 2013: selected articles from the eleventh asia pacific bioinformatics conference: bioinformatics. the full contents of the supplement are available online at http://www.biomedcentral.com/bmcbioinformatics/supplements/14/s 
our model consists of three dior-trinucleotide profiles identified through principle component analysis, and the predicted nucleosome occupancy flanking the poly sites. when we applied one model from one species to predict poly sites from other species, the sensitivity scores correlate with phylogenetic distances. a four-feature model geared towards small motifs was sufficient to accurately learn and predict poly sites across eukaryotes. we validated our model using two machine learning methods: logistic regression and linear discriminant analysis. we propose a poly site model that employs minimal features to capture the essence of poly sites, and yet, produces better prediction accuracy across diverse species. results show that models achieve 85-92% sensitivity and 85-96% specificity in seven animals and plants. we surveyed nine poly site prediction methods published between and all methods exploit the skewed nucleotide profile across the poly sites, and the highly conserved poly signal as the primary features for recognition.
 polyadenylation is present in all three domains of life, making it the most conserved post-transcriptional process compared with splicing and 5'-capping. even though most mammalian poly sites contain a highly conserved hexanucleotide in the upstream region and a far less conserved u/gu-rich sequence in the downstream region, there are many exceptions. furthermore, poly sites in other species, such as plants and invertebrates, exhibit high deviation from this genomic structure, making the construction of a general poly site recognition model challenging. we surveyed nine poly site prediction methods published between and all methods exploit the skewed nucleotide profile across the poly sites, and the highly conserved poly signal as the primary features for recognition. these methods typically use a large number of features, which increases the dimensionality of the models to crippling degrees, and typically are not validated against many kinds of genomes.  we propose a poly site model that employs minimal features to capture the essence of poly sites, and yet, produces better prediction accuracy across diverse species. our model consists of three dior-trinucleotide profiles identified through principle component analysis, and the predicted nucleosome occupancy flanking the poly sites. we validated our model using two machine learning methods: logistic regression and linear discriminant analysis. results show that models achieve 85-92% sensitivity and 85-96% specificity in seven animals and plants. when we applied one model from one species to predict poly sites from other species, the sensitivity scores correlate with phylogenetic distances.  a four-feature model geared towards small motifs was sufficient to accurately learn and predict poly sites across eukaryotes. 21- january the eleventh asia pacific bioinformatics conference vancouver, canada 