A web-based bioinformatics interface applied to the GENOSOJA project: databases and pipelines

RESEARCH ARTICLE

A web-based bioinformatics interface applied to the GENOSOJA project: databases and pipelines

Leandro Costa do Nascimento^I; Gustavo Gilson Lacerda Costa^I; Eliseu Binneck^II; Gonçalo Amarante Guimarães Pereira^I; Marcelo Falsarella Carazzolle^I,III

^ILaboratório de Genômica e Expressão, Departamento de Genética, Evolução e Bioagentes, Instituto de Biologia, Universidade Estadual de Campinas, Campinas, SP, Brazil

^IIEmpresa Brasileira de Pesquisa Agropecuária, Londrina, PR, Brazil

^IIICentro Nacional de Processamento de Alto Desempenho em São Paulo, Universidade Estadual de Campinas, Campinas, SP, Brazil

^{Send correspondence to}

ABSTRACT

The Genosoja consortium is an initiative to integrate different omics research approaches carried out in Brazil. Basically, the aim of the project is to improve the plant by identifying genes involved in responses against stresses that affect domestic production, like drought stress and Asian Rust fungal disease. To do so, the project generated several types of sequence data using different methodologies, most of them sequenced by next generation sequencers. The initial stage of the project is highly dependent on bioinformatics analysis, providing suitable tools and integrated databases. In this work, we describe the main features of the Genosoja web database, including the pipelines to analyze some kinds of data (ESTs, SuperSAGE, microRNAs, subtractive cDNA libraries), as well as web interfaces to access information about soybean gene annotation and expression.

Key words: bioinformatics, database, gene expression, soybean, Genosoja.

Introduction

Soybean is a legume of great economic importance in the international market, with a world production of almost 260 million tons for the 2009/2010 harvest. Brazil appears as the world's second largest producer, with a share of about 25%, and the crop is responsible for 10% of the country's total exports.

In recent years, the Brazilian soybean crop has been constantly threatened by climatic constraints (especially long periods of drought) and some attacks by pathogens such as Asian Rust disease (Yorinori et al., 2005), resulting in millions in losses for producers. Such constraints increase the importance of the breeding programs, established to discover new techniques for planting and prevention, increase production and lower costs. High-throughput sequencing technologies (like 454 pyrosequencing, Illumina/Solexa and ABI/SOLiD) represent significant advantages in these areas by producing millions of reads that can be used to measure levels of gene expression, allowing the identification of new genes or novel splice variants. On the other hand, it is necessary to intensify efforts by bioinformatics groups to develop new pipelines and data integration.

To this end, the Brazilian Soybean Genome Consortium (Genosoja) was established in 2007 with the goal of integrating several institutions currently working with soybean genomics in Brazil. The project promotes the search for solutions regarding possible treats, and to improve the soybean production process, emphasizing stresses that affect domestic production, like the occurrence of droughts, and pathogen attacks such as Asian Rust disease.

Among the main objectives of the Genosoja consortium is the creation of a relational database, integrating the results achieved by different methodologies and research groups working in the project. Despite the existence of other integrated soybean databases, such as SoyXpress (Cheng and Stromvik, 2008) and the "Soybean full-length cDNA database" (Umezava et al., 2008), as well as free pipelines available for data integration, such as Distributed Annotation System (DAS) (Dowell et al., 2001; Jenkinson et al., 2008) none of them integrates data from multiple experiments or provides transcriptome data from high throughput sequencing technologies like the database described in this work.

In light of this context, we created a soybean database connecting public soybean data (like ESTs and genomic sequences) and project data (like SuperSAGE tags, microRNAs and subtractive cDNA libraries). This database offers search tools for users, including keyword searches, statistical comparisons, automatic annotation (using some protein databases such as NR, Uniref, KEGG and Pfam), gene ontology classification and gene expression profiles under several conditions. Moreover, searches by sequence homology are possible using the local BLAST. All data are stored in a Fedora Linux machine, running the MySQL database server. The web interfaces (http://www.lge.ibi.unicamp.br/soybean) are based on a combination of CGI scripts using Perl language (including BioPerl module) and the Apache Web Server. As soon as the private data are published, the database will be freely available.

Methods, Results and Discussion

Public soybean data

In order to construct the Genosoja database we first collected all soybean data available at public biology databases. The genome of the cultivar Williams 82 and their predicted genes (66,153 sequences) were downloaded from the Phytozome (Schmutz et al., 2010). One full-length cDNA library from the Japanese cultivar Nourin2 was downloaded from the "Soybean full-length cDNA database". From NCBI (National Center for Biotechnology Information) we obtained 1,276,813 EST sequences (sequenced by SANGER and pyrosequencing technologies) and their equivalent GenBank files. All sequences were renamed in accordance to libraries, tissues and cultivars. This information was extracted from the GenBank files using homemade PERL scripts (Supplementary Material Figure S1). The bdtrimmer software (Baudet and Dias, 2005) was used to exclude ribosomal, vector, low quality and short (less than 100 bp) sequences. The EST assembly process was divided into two steps: (1) the ESTs were mapped into the soybean genome using the BLASTn algorithm (Altschul et al., 1997) (e-value cutoff of 1e^-10) and (2) all reads that aligned in same region of the reference were assembled together using the CAP3 program (Huang and Madan, 1999). The final result consists of 60,747 unigenes (30,809 contigs and 29,938 singlets). The effort to obtain the unigenes from assembled ESTs was important to increase the databases with information on untranslated regions (UTR), alternative splicing variants and gene expression profiling.

The Autofact program (Koski et al., 2005) was used to perform an automatic annotation of the predicted genes and the assembled unigenes. The main contribution of Autofact is the capacity to resume the annotation based on sequence similarity searches in several databases. For this, we used the BLASTx procedure (e-value cutoff of 1e^-5)to align the genes against certain protein databases, including: non-redundant (NR) database of NCBI, swissprot -databases containing only manually curated proteins (Suzek et al., 2007), uniref90 and uniref100 -databases containing clustered sets of proteins from UniProt, Pfam -a database of protein families (Bateman et al., 2002) and KEGG -a database of metabolic pathways (Kanehisa and Goto, 2000). The Autofact pipeline assigned function to 85% of the protein dataset. Figure 1 shows the complete pipeline of the public soybean data analysis.

Using the description of the origin of the ESTs (tissues and conditions), normalization procedures and statistical data analysis (Audic and Claverie, 1997), it was possible to infer differential gene expression among the assembled unigenes. This approach, called Electronic Northern, allows the users to compare gene expression profiles between two or more libraries and the results are available through a web interface (Figure 2).

Finally, the users can perform keyword and BLAST searches directly from the EST reads using the Gene Projects software (Carazzolle et al., 2007). This software also allows the user to perform assembly and annotation in these reads in an effort to improve unigene assembly. After generating a login/password it is possible to work on specific projects which users can develop and organize thematically by adding sequences to the assembly. After the assembly it is possible to view and to edit the results, improving the quality of the contigs.

Solexa SuperSAGE data

The Genosoja project generated three libraries using SuperSAGE methodology and these were sequenced by Illumina/Solexa technology. One library was constructed exploring gene expression in plants (Brazilian cultivar PI561356 -resistant) infected by the fungus Phakopsora pachyrhizi (Asian Rust disease) and two samples of plants (Brazilians cultivars: BR 16 -susceptible and Embrapa 48 -resistant) submitted to drought stress -for descriptions see Soares-Cavalcanti et al. (2012, this issue) and Wanderley et al. (2012, this issue). In total, the SuperSAGE approach generated 4,373,053 tags with 26 bp each.

Initially the tags of each sample were grouped in unique sequences. The unique sequences that presented low read counts (read count < 2) were discarded from the list. The Audic-Claverie statistic (Audic and Claverie, 1997) with a 95% confidence level (cutoff of 0.05) was used to identify tags as up-regulated (more expressed in the treated library) or down-regulated (more expressed in the control library).

In order to connect the unique tag with a gene sequence, the SOAP2 aligner program (Li et al., 2009) was used to align the unique tags with three databases (shown previously): (i) assembled unigenes (60,747 sequences), (ii) predicted genes (66,153 sequences) and (iii) the soybean genome. The program has been configured to allow for up to 2 mismatches in the alignments (SNPs can generate mismatches in the alignment, especially in this case because the assembled unigenes are generated by ESTs from different cultivars) and return all optimal alignments. In cases where more than one optimal alignment exists we decided to use the results from assembled unigenes, because they contain the UTR regions (a large part of the Super-SAGE tags are in the 3' UTR region), followed by the alignment with the predicted genes and, in the last case, the genome alignment was considered. Figure 3 shows the pipeline used in the SuperSAGE analysis.

For the sample submitted under drought stress (BR 16 and Embrapa 48 cultivars) we mapped 84.3% of the unique tags (Table 1), whereby the remaining 15.7% could represent new soybean genes specific to Brazilian cultivars. Similar results were obtained from a sample infected by Asian Rust, which yielded a total of 21,338 unique tags (20.42%) (Table 1) that did not align with any soybean database. In this case the unique tags from fungi genes may have contributed to increase this value. We tried to map these tags against Phakopsora pachyrhizi databases, but we did not find any fungus genes, probably due to the limited amount of fungus data available in the literature. The genome of this fungus, for example, has an estimated size of 500 Mb, but there is only 50 Mb available at NCBI.

We constructed a web interface for SuperSAGE analysis (Figure 4A). This interface shows, for each tag, the count number in both libraries (control and treated), the correspondent gene and its annotation (NR and Autofact result), as well as the position and the number of mismatches in the alignment. The user can filter the results using a keyword or gene name.

 

Solexa cDNA subtractive libraries data

Twenty-two cDNA subtractive libraries from different cultivars were sequenced in the Genosoja context, using many treatments with different time courses (Table 2) (Rodrigues et al., 2012, this issue). The reads were generated by Illumina/Solexa technology with read lengths of 45 or 76 bp, depending on the library.

In order to identify the genes in these libraries, the reads were mapped into soybean genes. First, we aligned the sequences against the unigenes using the SOAP2 aligner configured to allow up to two mismatches, discarding fragments with "Ns", and returning all optimal alignments. The sequences that did not align with unigenes were aligned against the predicted genes with the same parameters. A web interface (Figure 4B) provides users with all genes identified in each library and enables searches by gene name and keywords (in annotation results).

Solexa microRNA data

The Genosoja project generated eight small RNA libraries from soybean -four of the plants with Asian Rust disease (Brazilians cultivar PI561356 -resistant and Embrapa 48 -susceptible) and four under drought stress (Brazilians cultivars BR 16 -susceptible and Embrapa 48 resistant) (Molina et al., 2012, this issue). These libraries were sequenced using Illumina/Solexa technology and for each library the reads size range from 19 to 24 bp (Table 3).

Initially, the reads were grouped into unique sequences and read frequencies computed. The unique sequences that presented low read counts (read count = 2) were discarded from the list, as they were possibly caused by sequencing errors. In order to perform differential expression analysis between libraries, both a normalization and statistical significance analysis were applied using DEGseq software (Wang et al., 2009) considering a confidence level of 95% (cutoff of 0.05). Table 4 presents the number of unique and differential sequences in each library. For the statistical significance analysis, the treated over control libraries were considered.

To identify microRNAs from the small RNAs dataset it is necessary to identify the pre-microRNA by alignment of small RNAs (unique sequences) into the soybean genome assembly, followed by secondary structure identification. This alignment was performed using SOAP2 configured to allow for exact alignments only. The upstream and downstream genomic sequences of the read alignment position, 300 bp each in size, were extracted from the genome using homemade PERL scripts (Supplementary Material Figure S2). These genomic regions were aligned against the reverse complement of its respective tag (rc-tag) using the Smith-Waterman (Smith and Waterman, 1981) algorithm with two gaps and four mismatches allowed. The resulting sequences were considered pre-microRNA candidates, and the secondary structure was manually curated, resulting in 256 microRNAs (Figure 5) (Kulcheski et al., 2011).

Finally, the microRNA target prediction was performed using the Smith-Waterman algorithm (3 mismatches allowed) to align the 256 microRNAs against the assembled unigenes (shown previously). We considered only alignments in the 5'-3' direction obtained by comparison of the unigenes with the NR database using BLASTx. This methodology was able to identify targets for 169 microRNAs, most of which (39%) presented one or two targets (Figure 5).

Conclusions

In this work we presented all the bioinformatics analysis and pipelines used in the Genosoja project. The webbased interface constructed and described herein represents an important tool to help in the discovery of genes and new drugs that will enable increased soybean productivity. This system's use of common references (genome, assembled unigenes and predicted genes) facilitates the incorporation of new data from other sequencing methodologies or experimental conditions. Moreover, the bioinformatic pipeline discussed herein can also be applied to any genomic project, regardless of the organism.

Acknowledgments

The authors would like to acknowledge all research by the Brazilian Soybean Genome Consortium (Genosoja) involved in the generation of data used to construct the database presented in this paper. Moreover, we thank CNPQ (Conselho Nacional de Desenvolvimento Científico e Tecnológico -Brazil) for financial support to this work.

Internet Resources

NCBI, http://www.ncbi.nlm.nih.gov/ (October 10, 2011). Phytozome, http://www.phytozome.net/soybean.php (October 10, 2011). Soybean full-length cDNA database, http://rsoy.psc.riken.jp/ (Oc

tober 10, 2011). SOAP2 aligner, http://soap.genomics.org.cn/ (October 10, 2011). Genosoja database, http://www.lge.ibi.unicamp.br/soybean/ (Oc

tober 10, 2011). SoyXpress database, http://soyxpress.agrenv.mcgill.ca (October 10, 2011).

License information: This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Supplementary Material

The following online material is available for this article:

Figure S1 -Perl script to extract information about sequences from GenBank files.

Figure S2 -Perl script to extract the upstream and downstream genomic sequences of the read alignment position.

This material is available as part of the online article from http://www.scielo.br/gmb

^{Click to enlarge}

^{Click to enlarge}

• Altschul SF, Madden TL, Schaffer AA, Zhang J, Zhang Z, Miller W and Lipman DJ (1997) Gapped BLAST and PSI-BLAST: A new generation of protein database search programs. Nucleic Acids Res 25:3389-3402.
• Audic S and Claverie JM (1997) The significance of digital gene expression profiles. Genome Res 7:986-995.
• Bateman A, Birney E, Cerruti L, Durbin R, Etwiller L, Griffiths-Jones S, Howe KL, Marshall M and Sonnhammer ELL (2002) The Pfam protein families database. Nucleic Acids Res30:276-280.
• Baudet C and Dias Z (2005) New EST trimming strategy in: Brazilian Symposium on Bioinformatics, 2005. Lect Notes Bioinf 3594:206-209.
• Carazzolle MF, Formighieri EF, Digiampietri LA, Araujo MRR, Costa GLL and Pereira GAG (2007) Gene projects: A genome web tool for ongoing mining and annotation applied to CitEST. Genet Mol Biol 30(suppl):1030-1036.
• Cheng KCK and Stromvik MV (2008) SoyXpress: A database for exploring the soybean transcriptome. BMC Genomics 9:e368.
• Dowell RD, Jorkest RM, Day A, Eddy SR and Stein L (2001) The distributed annotation system. BMC Bioinformatics 2:e7.
• Huang X and Madan A (1999) CAP3: A DNA sequence assembly program. Genome Res 9:868-877.
• Jenkinson AM, Albrecht M, Birney E, Blankenburg H, Down T, Finn RD, Hermjakib H, Hubbard TJP, Jimenez RC, Jones P, et al. (2008) Integrating biological data -The Distributed Annotation System. BMC Bioinformatics 9(Suppl 8):S3.
• Kanehisa M and Goto S (2000) KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Res28:27-30.
• Koski LB, Gray LW, Lang BF and Burger G (2005) AutoFACT: An automatic functional annotation and classification tool. BMC Bioinformatics 6:e151.
• Kulcheski FR, Oliveira LFV, Molina LG, Almerao MP, Rodrigues FA, Marcolino J, Barbosa JF, Stolf-Moreira R, Nepomuceno AL, Marcelino-Guimaraes FC, et al. (2011) Identification of novel soybean microRNAs involved in abiotic and biotic stress. BMC Genomics 12:e307.
• Li R, Yu C, Li Y, Lam T-W, Yiu S-M, Kristiansen K and Wang J (2009) SOAP2: An improved ultrafast tool for short read alignment. Bioinformatics 25:1966-1967.
• Molina L, Cordenonsi G, Loss G, Oliveira LFV, Carvalho K, Kulcheski F and Margis R (2012) Metatranscriptomic analysis of small RNAs present in soybean deep sequencing libraries. Genet Mol Biol 35(suppl 1):292-303.
• Rodrigues FA, Marcolino J, Carvalho JFC, Nascimento LC, Neumaier N, Farias JRB, Carazzolle MF, Marcelino FC and Nepomuceno AL (2012) Using subtractive libraries to prospect differentially expressed genes in soybean plants submitted to water deficit. Genet Mol Biol 35(suppl 1):304-314.
• Schmutz J, Cannon SB, Schlueter J, Ma J, Mitros T, Nelson W, Hyten DL, Song Q, Thelen JJ, Cheng J, et al. (2010) Genome sequence of the paleopolyploid soybean. Nature 463:178-183.
• Smith TF and Waterman MS (1981) Identification of common molecular subsequences. J Mol Biol 147:195-197.
• Soares-Cavalcanti NM, Belarmino LC, Kido EA, Pandolfi V, Marcelino-Guimaraes FC, Rodrigues F, Pereira GAG and Benko-Iseppon AM (2012) Overall picture of expressed Heat Shock Factors in Glycine max, Lotus japonicus and Medicago truncatula. Genet Mol Biol 35(suppl 1): 247-259.
• Suzek BE, Huang H, McGarvey P, Mazumber R and Wu CH (2007). Uniref: Comprehensive and non-redundant UniProt reference clusters. Bioinformatics 23:1282-1288.
• Umezawa T, Sakurai T, Totoki Y, Toyoda A, Seki M, Ishiwata A, Akiyama K, Kurotani A, Yoshida T, Mochida K, et al. (2008) Sequencing and analysis of approximately 40,000 soybean cDNA clones from a full-length-enriched cDNA library. DNA Res 15:333-346.
• Wang L, Feng Z, Wang X, Wang X and Zhang X (2009) DEGseq: An R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics 26:136-138.
• Wanderley-Nogueira AC, Belarmino LC, Soares-Cavalcanti NM, Bezerra-Neto JP, Kido EA, Pandolfi V, Abdelnoor RV, Binneck E, Carazzole MF and Benko-Iseppon AM (2012) An overall evaluation of the Resistance (R) and Pathogenesis Related (PR) superfamilies in soybean, as compared with Medicago and Arabidopsis Genet Mol Biol 35(suppl 1):260-271.
• Yorinori JT, Paiva WM, Frederick RD, Costamilan LM and Bertagnolli PF (2005) Epidemics of soybean rust (Phakopsora pachyrhizi) in Brazil and Paraguay from 2001 to 2003. Plant Disease 89:675-677.

Publication Dates

• Publication in this collection
  01 June 2012
• Date of issue
  2012