Joly et al. BMC Genomics 2010, 11:422 http://www.biomedcentral.com/1471-2164/11/422 digested with NlaIII and the band containing the ditags was extracted from 12% acrylamide gel. The purified ditags were self-ligated to form concatemers of 500- 1800 bp isolated by agarose gel. The resulting DNA fragments were ligated into the SphI site of pUC19 and cloned into UltraMAX DH5aFT E. coli cells (Invitrogen, Carlsbad, CA, USA). White colonies were screened by PCR to select long inserts for automated sequencing as previously described for cDNA libraries. Sequence files were analyzed using the SAGEparser program [83]. Tags corresponding to linker sequences were discarded and duplicate concatemers were counted only once. To identify the transcripts, the sequences of 15 bp SAGE tags (NlaIII site CATG plus adjacent 11 bp tags) were matched with a collection of unassembled poplar and M. larici-populina ESTs with polyA tail using a customized Python script. To avoid the possibi- lity of sequencing errors in the EST database, the matches that were identified only once among the EST database were not considered. Positive selection analyses Homologous gene groups (HGGs) were identified from reciprocal TBLASTX searches (E-value < 1e-30) between libraries and P. graminis f. sp. tritici gene models fol- lowed by graph clustering analysis using a TCL imple- mentation of the Deep-First Search algorithm. Each HGG included at least two sequences from Melampsora unisequences and a minimum of three total sequences. HGGs were either removed or divided when overlapping regions where too short or when similarity was not found throughout the majority of the unisequence coding sequences. Furthermore, sequences with gaps across the aligned coding sequences were removed in order to mini- mize the impact of the pitfalls of positive selection ana- lyses, such as gap-induced misalignments and relaxed selection in pseudogenes. The resulting 369 HGGs were then submitted to positive selection analyses using a suite of program grouped into a single Python script. The pro- tein sequences in each HGG were first aligned with Clus- talW [84] and the corresponding coding DNA sequences were automatically extracted. A Neighbor-joining phe- netic tree based on distance matrix between nucleotidic sequences was then reconstructed for each HGG and used as starting tree for Bayesian inference and Markov Chain Monte Carlo simulations (B/MCMC) (only possi- ble with HGG of 4 sequences and more; the Neighbor- joining tree reconstructed with PAUP* was used for HGG of 3 sequences). Prior to the B/MCMC, the models for nucleotide substitutions were selected using the hier- archical likelihood ratio test (hLRT) implemented in the Modeltest 3.7 program. Adaptive evolution was first esti- mated by pairwise calculation of the rates of nonsynon- ymous nucleotide substitutions per nonsynonymous site Page 13 of 16 (dN) and the rates of synonymous nucleotide substitu- tions per synonymous site (dS) between all members of an HGG. Additionally, as adaptive evolution is likely to act on a small subset of amino acid residues and hence averages of substitution rates across the gene may not strictly indicate positive selection, HGGs were scanned for adaptive evolution using codon-based substitution models that allow ω to vary among sites, with the para- meters of the model estimated using maximum likeli- hood. These analyses were conducted using the codeml application from the PAML package version 4 [67]. Baye- sian inference of phylogeny aimed at estimating the pos- terior probabilities and branch length of phylogenetic trees as starting values for codeml maximum likelihood iteration under codon model M0 to get fixed branch lengths. The resulting ω pairwise calculations are shown in Table 3 for all HGGs with at least one pairwise calcu- lation showing ω > 1.2. We also contrasted the codon substitution models M1A (neutral), M2A (selection), M7 (beta), M8 (beta and ω) and M8A (beta and ω = 1; [85]). The model M1A assumes two site classes in proportions p0 and p1 = 1-p0 with 0 < ω 0 < 1 (conserved) and ω 1 = 1 (neutral). M2A adds an additional class of site with ω 2 as a free parameter, allowing for sites with ω 2 > 1 with proportion p2. Model M7 uses a beta distribution of sites within the interval 0 < ω < 1. M8 adds an extra class of sites to the M7 model, allowing for positively selected sites with ω > 1, while this extra dN/dS category is restricted to one in model M8A [85]. From these models, statistical significance was tested using likelihood ratio tests by comparing the models M1A, M7 and M8A with the alternative M2A, M8 and M8 models, respec- tively. Models M2A and M8 are tests of positive selection among codon sites and were implemented with at least three different starting ω values (0.2, 1.0 and 2.0). Twice the difference in log likelihood ratio between and alternative models was compared with a c2 distribution with two degrees of freedom. The likelihood ratio tests assess whether the alternative models fits the data better than the model and is known to be conservative in simulation tests. For HGGs that tested positive using the ML method, posterior Bayesian probabilities of site classes were inferred for each amino acid site by using the Bayes empirical Bayes method [86]. Additional material Additional file 1: PFAM domains found in the putative secreted members of the Melampsora unisequence dataset. Domains represented in more than five unisequences are shown. PFAM hits were considered significant when E-value ≤ 1e-5. Additional file 2: Top five most abundantly represented secretome members in each Melampsora library. BLASTX hits were considered significant when E-value ≤ 1e-4.