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 .
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 .
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 .
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 ) and ( 2 ) all reads that aligned in same region of the reference were assembled together using the cap3 program ( huang and madan , 1999 ) .
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 ) 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 ) .
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 .
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 .
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 .
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
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 .
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 .
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 ,
first , we aligned the sequences against the unigenes using the soap2 aligner configured to allow up to two mismatches , discarding fragments with 
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 ) .
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 .
these libraries were sequenced using illumina / solexa technology and for each library the reads size range from 19 to 24 bp ( table 3 ) .
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
table 4 presents the number of unique and differential sequences in each library . for the statistical significance analysis ,
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 .
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 ) .
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 .
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 ) and ( 2 ) all reads that aligned in same region of the reference were assembled together using the cap3 program ( huang and madan , 1999 ) .
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 ) 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 ) .
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
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 .
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 .
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
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 .
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 .
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 ,
first , we aligned the sequences against the unigenes using the soap2 aligner configured to allow up to two mismatches , discarding fragments with
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 ) .
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 ) .
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
table 4 presents the number of unique and differential sequences in each library . for the statistical significance analysis ,
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 .
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 ) .
in this work we presented all the bioinformatics analysis and pipelines used in the genosoja project . the web - based 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 .
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 .
 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 . 
Introduction
Methods, Results and Discussion
Public soybean data
Solexa SuperSAGE data
Solexa cDNA subtractive libraries data
Solexa microRNA data
Conclusions
Supplementary Material