Wang et al. BMC Genomics 2012, 13:529 http://www.biomedcentral.com/1471-2164/13/529 empty reads and low quality sequences (reads with un- known sequences 'N') [30]. Next, the reads were ran- domly clipped into 21 bp K-mers for assembly using SOAPdenovo software because the 21-mer provided the best result for transcriptome assembly. Small K-mers made the graph very complex; while large K-mers have poor overlap in low sequencing depth regions [31]. For alleles, the nucleotides with the highest frequency were selected. The resultant contigs were joined into scaffolds based on the mate pairs information and the scaffolds were clustered using TGI Clustering tools [47]. Assembled genes were used for subsequent analyses and are referred to as “distinct sequences”. Data deposit The data sets of Illumina sequencing are available at the NCBI Short Read Archive (SRA) database with the acces- sion number: SRR062575. The assembled sequences were deposited in the NCBI Transcriptome Shotgun Assembly (TSA) database under the accession number of HP777244 to HP823074 and can be searched using the GeneID listed in Additional file 2, Supporting information. Functional annotation and gene expression analysis Distinct sequences were annotated by Blast search against the NCBI nr database with a cut-off E-value of 10-5. GO annotation was analyzed by Blast2GO software [48]. The GO terms were retrieved from Blast hits with an E-value threshold of 10-5. Comparisons of the distri- bution of GO terms among the Asia II, MEAM1 and MED transcriptomes were done using the Web Gene Ontology Annotation Plot (WEGO) [49]. Pathway anno- tation was performed using Blastall software against the KEGG database. Based on the number of reads for a gene, gene expression levels can be estimated from Illu- mina sequencing with great accuracy [50,51]. To esti- mate the level of gene expression, the number of reads mapped to each distinct sequence was extracted. Since read mapping is sensitive to the size of the target refer- ence sequence and sequencing amount, we adjusted the raw read count by the total number of reads mapped and the length of the gene by calculating Reads Per Kilo- base per Million mapped reads (RPKM) [50]. Analysis of protein families To reveal the functional differences among Asia II 3, MEAM1 and MED transcriptomes, we analyzed their protein families. Using Blastx (E-value <10-5), the trans- lated region of each gene was identified by aligning the sequence to the Swissprot database. The longest trans- lated protein sequence for each gene was then extracted and sequences less than 200 bp were removed. To iden- tify protein families among the three transcriptomes, an all-against-all Blastp was performed for all the translated Page 10 of 13 genes from the three transcriptomes. Blastp results were further analyzed by a Markov Cluster Algorithm (MCL) with an inflation factor of 1.6. The protein families belonging to all the three transcriptomes were referred to as “Core protein families”. Based on the GO annota- tion of Asia II 3, MEAM1 and MED transcriptomes, we calculated the total number of genes under each GO term in the “Core protein families” and the three tran- scriptomes. For each GO term, its enrichment in the “Core protein families” was measured using the hyper- geometric test with an cut-off p value of 10-5 [33]. Identification of orthologous genes and prediction of coding and untranslated regions The orthologous genes between Asia II 3 and MEAM1 and those between Asia II 3 and MED were identified respectively according to the previous description using MegaBLAST [30]. Briefly, pairs of sequences that were reciprocally a best hit and with a minimum length of 200 bp were retained as putative orthologs. To remove potential paralogs, only pairs of sequences unambigu- ously mapped to the same protein in Swissprot database with an E-value <1×10-5 were selected. The CDS of the orthologous genes were determined by BLASTx the Swissprot database with an E value<1×10-5. The start codon was determined by examination of the in-frame ATG codon of the aligned reference protein. The stop codon position was determined by examination of in- frame TAA, TAG and TGA motifs present within 30 bp of the stop codon of the reference protein. The 5’UTR and 3’UTR regions were defined based on the position of start codon, stop codon and predicted CDS. To pre- vent false positive results, only UTR pairs with an E- value<1×10-30 were selected for further analyses. CDS containing unexpected stop codon(s) and shorter than 150 bp were removed. Sequence divergence analyses and estimation of substitution rates The 5’UTR, CDS and 3’UTR regions were separately extracted from each pair of orthologs. The CDS and UTR regions were aligned separately to each other with a MegaBlast algorithm and checked manually for errors. Only the homologous regions of each gene pair were extracted for sequence comparison. Sequence divergence between the homologous regions of each gene pair was calculated by dividing the number of substitutions by the number of base pairs compared. The average diver- gence between transcriptomes was determined by divid- ing the total number of substitutions by the total number of base pairs compared. The sequence diver- gence at nondegenerate (nd), fourfold degenerate (4d), CpG and non-CpG regions was determined respectively according to the previous descriptions [39]. The ratio of