BACKGROUND
comparative studies of gene expression have been used to understand the regulation of a wide range of biological processes.
with the development of next generation sequencing  technologies, and in particular the use of illumina sequencing platforms, reliable genome wide comparison of gene expression between different biological conditions has become possible .
moreover, a growing number of available genome and transcriptome sequences  now provides the opportunity to compare gene expression not only in well-established, but also in emerging model systems.
especially, the comparison of gene expression between both closely  and distantly related species  has great potential to help understand phenotypic divergence and species adaptations at a mechanistic level  <cit> .
experiments to study differential gene expression using ngs technologies  are based on a sequencing library generated from reverse transcribed messenger rna  that is extracted from the tissue and conditions of interest.
illumina sequencing, for example, results in the generation of millions of short reads ranging from 36 bp to 150 bp  <cit> .
the first step of the bioinformatics analysis is to align these reads to a reference that represents all transcripts that should be quantified .
this reference can be a whole genome sequence with annotated gene models or a transcriptome.
the latter can either be generated by a de novo assembly of the rna-seq reads  <cit>  or it could be extracted from an annotated genome.
the next step is to determine the number of reads that are aligned to a gene model or transcript.
depending on the type of reference used  various different methods have been established  <cit> .
finally, the number of reads assigned to a given gene model or transcript is compared between different conditions to identify differentially expressed genes.
the steps outlined above for a general rna-seq experiment are suitable to compare gene expression levels between different conditions, stages or tissues of the same species.
however, comparison of gene expression between different species or populations of the same species needs to account for differences in gene sequences.
in this case, reads should be mapped to species-specific references for which the expression level of a gene in one of the species is compared to the expression level of its ortholog in the other species.
most importantly, this requires sets of orthologous genes reliably identified in all references.
since genomes or transcriptomes are usually generated by different research groups for different applications and using different pipelines for assembly and annotation, annotated references for inter-specific gene expression studies are often not comparable.
for instance, orthologous genes might be missing from one or more of the references as result of natural variation or technical problems like incomplete assemblies or too many sequencing errors, which hampers unequivocal identification of orthologous genes.
additionally, it is common practice to filter out genes that are incomplete or lack synteny in relation to a model reference from new gene model predictions  <cit> .
even though there are many tools available to perform genome annotation, a general standard does not exist.
therefore, the final gene set generated by each genome project will have genes missing as a result of methodological problems and filtering criteria, and this can directly influence the result of the differential gene expression analysis  <cit> .
additionally, even if most one-to-one orthologs have been successfully identified in different references, these gene models may vary in length for various reasons: first, the genes could naturally differ in length among species.
second, as a consequence of differences in the sequence or assembly quality of the reference genomes , orthologous gene models might be truncated in one or more of the references.
to our knowledge, a comprehensive evaluation of methods that could be applied to account for inter-specific gene length differences has not been performed yet.
a plethora of statistical approaches have been developed to determine whether differences in the number of reads are due to technical variation or due to real biological differences in gene expression.
detailed evaluation and comparison of these methods concluded that the most accurate statistical validation of differential gene expression is reached when statistical models are used that directly take the number of aligned reads into account .
these methods include standard and generalized poisson and negative binomial distributions to model count-based expression data  <cit>  as implemented in deseq  <cit> , deseq <dig>  <cit> , edger  <cit>  or degps  <cit> .
also the differential expression analysis based on moderated t-statistics as implemented in the limma package  <cit>  using log-transformed count per million values originating from normalization with voom  <cit>   has been shown to perform extremely well  <cit> .
while all of these methods account for most technical biases and control well false positive rates, none of these methods is specifically designed to account for gene length differences as they occur in inter-specific expression studies.
one potential solution could be the application of the normalization method reads per kilobase per million mapped reads  as it accounts for length differences in gene models  <cit> .
however, it has been shown that even after correcting for length differences, a longer transcript is more likely to appear as differentially expressed if rpkm values are used to assess the statistical significance .
rpkm normalization is still widely used to compare gene expression levels of different genes within a species, but to our knowledge it has not been tested if this method efficiently normalizes length differences when comparing gene expression in different species.
here we show that the published genomes of three closely related drosophila species, d. melanogaster, d. simulans and d. mauritiana have qualitative limitations as references for comparative gene expression studies.
this is mainly due to the fact that many genes cannot be properly compared because orthologous genes are missing in the annotation of at least one of the genomes.
even after a direct re-annotation of the three genomes using the same annotation pipeline many orthologous gene models exhibit significant length differences.
taking advantage of these inter-specific gene length differences in the published and the directly re-annotated references, we benchmarked four statistical frameworks  for their ability to reduce the number of potentially false positives.
we demonstrate that none of these methods sufficiently accounts for the observed differences in gene length.
therefore, we propose that the length normalization should be performed prior to read mapping during the generation of the mapping references.
we report a straightforward re-annotation method that relies on a reciprocal re-annotation of orthologous gene models in two or more species.
this approach allows the comparison of nearly all genes that have been annotated in d. melanogaster in all three species.
additionally, we find that the use of the new gene sets as mapping references results in a more robust estimation of transcript abundance and a more reliable comparison of gene expression levels between species.
we propose that the generation and annotation of new genome resources or the re-annotation of existing genomes will be powerful tools to establish gene expression profiling in many emerging model systems.
RESULTS
analysis of published genome annotations reveals a reduced number of comparable gene models for differential gene expression studies between species
we first assessed the completeness and comparability of the published gene sets for the three closely related species d. melanogaster, d. simulans and d. mauritiana.
at the time of our analysis, the annotation of the d. melanogaster genome  - one of the best curated metazoan genomes available at flybase  - included  <dig>  unique protein coding genes.
the most recent annotations for d. simulans  <cit>  and d. mauritiana  <cit>  were generated using the d. melanogaster gene set as reference .
both gene sets contain a large fraction of the  <dig>  d. melanogaster genes .
however, orthologs of almost  <dig> d. melanogaster genes are not included in each of the final gene sets either because the authors applied various filtering steps to exclude incomplete orthologous sequence with respect to the d. melanogaster gene  or because the genes are not present in one of the species.
since these filtering steps are influenced by the quality of each of the assembled genome and the scientific question of each research group, the missing genes in both annotations are not the same.
only  <dig> genes  can be identified unequivocally as orthologs in all three annotations .
among the genes missing in at least one annotation, we found some important and well-studied developmental genes including the hox genes abdominal b , ultrabithorax  or antennapedia , the head and brain patterning gene orthodenticle  and the segment polarity gene hedgehog  .table  <dig> number of genes obtained by each annotation method
d. melanogaster
d. simulans
d. mauritiana
the last column contains the number of genes for which 1: <dig> orthologs were identified in the three species.
“after filtering” indicates the remaining common genes after filtering out genes with length difference larger than  <dig> bp.
percentages in brackets are always given in relation to the total number of gene models in d. melanogaster
next we assessed the comparability of the three reference genome annotations in terms of gene length, since length differences larger than the length of the rna-seq reads are likely to introduce a bias during mapping and the subsequent differential expression analysis.
if we consider 50 bp single-end reads, which have been shown to be long enough to produce accurate results when measuring differential gene expression , genes that have a length difference larger than 49 bp among the annotations of the three drosophila species are likely to bias a subsequent differential gene expression analysis.
a pairwise comparison of annotated gene length for the  <dig> genes present in all three drosophila species  shows that in the published annotations, the gene length differences are larger than 49 bp in  <dig>  %  of the orthologous genes between d. mauritiana and d. simulans,  <dig>  %  between d. melanogaster and d. simulans and  <dig>  %  between d. melanogaster and d. mauritiana .
if these length differences are not accounted for, these genes could result in false positives in a differential gene expression analysis.fig.
<dig> pairwise length difference between orthologous genes.
bars indicate the number of genes with that difference in length  for each pair of species.
green shades indicate differences lower than 50 bp while orange to red indicate larger differences.
the comparison is showed for  the published annotations, b the direct re-annotation of the published genomes and c the reciprocal re-annotation of the published genomes
direct re-annotation of published genomes
next we asked whether a direct re-annotation of the d. simulans and d. mauritiana genomes individually using the same d. melanogaster gene set as reference and the same annotation pipeline allows the comparison of more genes in an inter-specific differential gene expression study.
we used the  <dig>  protein sequences of d. melanogaster  as reference to re-annotate the published genomes of d. simulans  <cit>  and d. mauritiana  <cit>  using the program exonerate  <cit> .
without applying any filtering, we find orthologs of  <dig>  d. melanogaster genes that are comparable among the two species .
next, we determined the length of directly re-annotated genes that are found in all three species.
this comparison shows an increase in number of genes with the same length between the three species after direct re-annotation .
however, a high number of orthologous genes have a length difference of more than 49 bp :  <dig>  between d. melanogaster and d. mauritiana,  <dig>  between d. melanogaster and d. simulans and  <dig>  between d. mauritiana and d. simulans.
these observed length differences could be due to real natural variation in coding sequence length between species or they could be technical artifacts, for example truncated gene sequences arising from sequencing or genome assembly errors.
in summary, although annotated genomes are available for the three closely related drosophila species d. melanogaster, d. simulans and d. mauritiana, their use as mapping references for inter-species differential gene expression analyses is limited due to missing orthologs and a potential bias because of different annotated gene lengths.
the use of the same reference gene set, annotation pipeline and the lack of filtering incomplete gene sequences results in an increase in the number of comparable genes in these three closely related species.
length difference in reference genes introduces biases in differential expression studies
since we find a high number of gene models with length differences > 49 bp in the published annotations and after the direct re-annotation , the three drosophila genomes are excellent models to test whether length differences larger than the read length do indeed influence the statistical analysis of differential gene expression.
we used the published d. melanogaster annotation and the newly generated direct re-annotations of d. simulans and d. mauritiana as mapping references for pairwise comparisons of gene expression between d. melanogaster and d. mauritiana and d. simulans and d. mauritiana using 50 bp single-end illumina rna-seq reads generated for these three species .
rna-seq reads generated from one species were mapped only to the gene set of that species.
using this experimental setup, we compared four different statistical frameworks, namely deseq <dig>  deseq <dig> with gene length correction  <cit> , limma with length correction based on rpkm  <cit>  and voom-limma  <cit>  including rpkm length correction.
for each method, we first report the number of differentially expressed genes for each of the two pairwise species comparisons.
next we evaluate the impact of the length differences between gene models on the fold-change in gene expression between species.
and eventually, we compare the results of each of the four methods to an independent qpcr experiment for a subset of genes.
deseq <dig> without length correction
first we performed the statistical analysis for differential gene expression with deseq <dig>  <cit>  using directly the read counts for each gene model.
for both pairwise comparisons using the published and the direct re-annotation as reference, we found that between  <dig>  and  <dig>  % of all comparable genes are significantly differentially expressed .
additionally, we found a very strong correlation between inter-specific length differences of the gene models  and the log2-fold change in gene expression .
the negative correlation means that genes that are longer in one species appear to be more up-regulated and vice versa.
the correlation can be explained by the mapping procedure: in orthologous genes with different length, more reads align to the ortholog that has the longer gene model .
this results in an artificially higher expression for this specific gene in the species with the longer gene model.
from this correlation we also see that most of those genes with length differences and a high log2-fold change are also significantly differentially expressed , showing that this method introduces a large number of false positives.fig.
<dig> length differences between orthologous genes introduce gene expression biases.
relation between length differences and the log2-fold change in the rna-seq experiment between d. mauritiana and d. simulans using the direct re-annotation of these species as mapping references.
dots represent genes with length difference > 49 bp in these annotations .
genes significantly differentially expressed in the presented analysis  are shown in red.
a negative log2-fold change indicates higher expression in d. mauritiana.
a positive length difference indicates that the ortholog of d. mauritiana is longer.
the p-value and rho of the spearman’s rank correlation are indicated on the lower right side of the plots.
a results of deseq <dig> without length correction.
b results of deseq <dig> applying length normalization factor matrix.
c results of applying rpkm normalization and limma to call differentially expressed genes.
d results of applying voom normalization followed by a length normalization matrix and limma to call differentially expressed genes
fig.
<dig> schematic representation of length bias in inter-species differential expression analysis and our reciprocal re-annotation strategy to correct it.
a length bias in the analysis of a non-differentially expressed gene.
coloured rectangles represent the part of the transcript which is included as reference for the rna-seq reads to map to, while unfilled rectangles are regions of the transcript which are omitted and to which rna-seq reads cannot be mapped.
red “n”s represent sequencing errors that prevent the complete annotation of a transcript.
mapped reads are shown as thin black lines and the number bellow indicates the total of reads mapped.
if one transcript is shorter in one of the references compared to its orthologs, for the same expression levels fewer reads will map to it.
this can result in false positives in the analysis of differential expression.
our strategy to correct this bias is to shorten the orthologs in the other references to match the length of the shorter sequence.
b pipeline of reciprocal transcriptome re-annotation method.
black numbers in white circles represent genome annotation steps using the “est2genome” command of exonerate  <cit> .
grey numbers in grey circles represent conversion of the resulting gff file into a new transcript set.
filled horizontal bars represent the annotated set of transcripts; non-filled horizontal bars at the start/end of the transcripts represent parts of the transcript that cannot be correctly annotated in one reference and are therefore eliminated from the transcript set.
the boxes with red frame indicate the transcript sets that will be used as reference for rna-seq read mapping .
step 1: the transcript set of the best annotated genomes  is used to annotate one of the other genomes  and generate a new transcript set for this species.
due to sequencing errors, some transcripts will be shorter.
step 2: the new transcript set form d. simulans is used to annotate the last genome .
the gene set generated contains shorter transcripts due to sequencing errors in d. mauritiana but also in d. simulans.
step 3: the transcript set from d. mauritiana is used to re-annotate the previously generated set from d. simulans to integrate the information from the d. mauritiana assembly.
step 4: the second transcript set from d. simulans is used to annotate the d. melanogaster set in order to integrate the information from d. simulans and d. mauritiana
in order to specifically test whether differences in the length of gene models indeed influence the differential expression analysis we chose seven genes that were shorter in the d. mauritiana published annotation compared to the published d. melanogaster annotation.
when we analysed the differential expression using deseq <dig> without any length correction, the expression of all seven genes were significantly different .
the log2-fold change value indicated that d. melanogaster had a significantly higher expression than d. mauritiana .
to validate the results obtained by the rna-seq experiment, we measured the relative expression of the seven genes in d. melanogaster and d. mauritiana using qpcr.
this analysis showed that the seven genes that had length differences in the species-specific annotations were not significantly differentially expressed .
as a control we chose another three genes that showed significant differential expression in the rna-seq data but had the same length in both species in the two annotation methods .
we found that piwi and alrm showed a significantly higher expression in d. mauritiana when using this alternative quantification method, confirming the results obtained by rna-seq .
moreover, nplp <dig> had higher expression in d. melanogaster again consistent with our rna-seq data, although this difference was not significant .table  <dig> differentially expressed genes and correlation between calculated log2-fold changes and length difference between orthologous genes.
results are shown for the four applied methods, the three studied annotation strategies and the two described pairwise species comparisons
d. mau vs d. mel
d. mau vs d. sim
d. mau vs d. mel
d. mau vs d. sim
d. mau vs d. mel
d. mau vs d. sim
d. mau vs d. mel
d. mau vs d. sim
a
d. mau vs d. mel
a
d. mau vs d. sim
d. mau vs d. mel
d. mau vs d. sim
d. mau vs d. mel
d. mau vs d. sim
d. mau vs d. mel
a
d. mau vs d. sim
d. mau vs d. mel
d. mau vs d. sim
d. mau vs d. mel
d. mau vs d. sim
a
d. mau vs d. mel
a
d. mau vs d. sim
a
d. mau vs d. mel
d. mau vs d. sim
afdr  <dig>
bspearman’s rank correlation is measured between log2fc and length difference of genes showing more than 49 bp length difference
published annotation: d. mau vs. d.mel:  <dig>  genes/d.mau vs. d.sim:  <dig> genes; direct annotation: d.mau vs. d.mel:  <dig> genes/d.mau vs. d.sim:  <dig> genes; reciprocal annotation: d.mau vs. d.mel:  <dig> genes/d.mau vs. d.sim:  <dig> genes
*p <  <dig> ; **p <  <dig> ; ***p <  <dig>
fig.
<dig>
qpcr results.
boxplot of normalized ct values  for each studied gene  boxplot is showed for ct values in d. melanogaster orer  and d. mauritiana tam <dig> .
or t-welch-test ; *p <  <dig> , **p <  <dig> ; ***p <  <dig> )
in summary, these results suggest a high level of potentially false positive candidates when methods based on direct read counts without the application of length correction are used with mapping references that exhibit differences in the length of orthologous genes.
deseq <dig> with length correction
next we benchmarked the use of deseq <dig>  <cit>  including a normalization factor matrix incorporating gene length to account for the length differences between orthologous genes.
using this approach for pairwise gene expression analyses, we found that only  <dig>  to  <dig>  % of the comparable genes were significantly differentially expressed.
even though the gene length was accounted for during the deseq <dig> analysis of differential gene expression, we still find a correlation between inter-specific gene length differences and log2-fold changes .
however, the significance of this correlation is greatly reduced in comparison to the deseq <dig> analysis without length correction , suggesting that the length correction incorporated in deseq <dig> helps to reduce the number of false positive candidate genes.
this finding is further supported by the comparison of the rna-seq results to the qpcr data.
after length correction only one  of the seven genes that are longer in d. melanogaster show significant differential expression .
limma with rpkm length correction
rpkm values are commonly calculated for rna-seq datasets to account for variation in library sizes and to correct for length differences between different genes within the same species  <cit> .
the moderated t-statistics incorporated in the limma r package  <cit>  can subsequently be used to assess differential gene expression.
it has not been tested if this approach also corrects properly for differences in the length of the same gene being compared between two species.
using this method, we found between  <dig>  % and to  <dig>  % of the comparable genes to be significantly differentially expressed .
our correlation analysis shows that the correction of a length bias using rpkm values still results in a clearly significant correlation between the gene length difference and the observed log2-fold change .
however, compared to the deseq <dig> analysis without length correction, the significance values are highly reduced , showing that the number of false positives is lower.
accordingly, six of the seven genes that we benchmarked with qpcr show no significant differential gene expression although they show clear length differences between d. melanogaster and d. mauritiana .
again cp <dig> is the only gene that appears as significantly differentially expressed also after correcting for length differences.
voom-limma with rpkm length correction
it has recently been shown that differential gene expression analysis with limma  <cit>  using normalized read counts from voom  <cit>  perform very well for rna-seq datasets  <cit> .
although this method is designed to work with direct read counts, in this case we tested it with an additional transcript length correction.
between  <dig> and  <dig>  % of the comparable genes are significantly differentially expressed .
after length correction  and normalization with voom, we found a significant correlation between gene length differences and log2-fold changes when the published annotations and the directly re-annotated reference gene sets were used however, this was slightly reduced compared to the rpkm-limma analysis, especially for the d. simulans and d. mauritiana comparison.
.
for the seven qpcr benchmarked genes that have clear length differences between d. melanogaster and d. mauritiana the use of the voom-limma method results in three significantly differentially expressed genes , suggesting a higher false positive rate.
length correction during the statistical analysis might be insufficient
the comprehensive comparison of four methods for differential gene expression analysis shows that the incorporation of a length correction drastically reduces the number of false positive candidate genes.
although the correlation between gene length differences and the observed log2-fold changes  is reduced in the three methods that account for gene length differences , none of them sufficiently corrects the length bias present in the two gene sets used as mapping references.
this is also supported by the qpcr validation of seven genes that exhibit clear length differences between the published d. melanogaster and d. mauritiana annotations .
in all three methods at least one gene was still artificially significantly differentially expressed.
this is most pronounced for the voom-limma method where three of the seven genes are significantly differentially expressed.
of the seven genes we analysed using qpcr, cp <dig> was in all cases identified as a false positive candidate.
in order to further characterise this gene, we visually inspected the distribution of mapped reads.
interestingly, the 3′-region is missing in the d. mauritiana ortholog of cp <dig>  and we found clearly more d. melanogaster reads that map to this 3′-part of the transcript than to the 5′-region .
the number of d. melanogaster and d. mauritiana reads mapped to the 5′-part of the transcript is comparable .
hence, a very likely explanation for the inefficient length correction of the three applied methods could be an unequal distribution of the mapped reads along the transcripts.
besides an insufficient length correction, the deseq <dig> method including a length matrix results in the lowest number of significantly differentially expressed genes, suggesting that the length correction applied here might be extremely conservative and could lead to a high rate of false negatives.
interestingly, in many pairwise comparisons we found more significantly differentially expressed genes with a higher expression in d. mauritiana compared to d. melanogaster and d. simulans , although d. mauritiana gene models are generally shorter than those of the other two species .
this finding suggests that the length correction applied here might reduce the power to detect differential expression for the already short d. mauritiana genes.
in summary, all three methods that include a length correction decrease the chance of identifying false positives.
the rpkm-voom-limma and rpkm-limma methods seem to give the best ratio of false positives and false negatives, while deseq <dig> including a length matrix is very conservative.
however, none of the length correction methods tested does efficiently account for all differences in gene length observed in the reference annotation of the three studied drosophila species.
the length bias is most obvious when the distribution of reads is not uniform across the transcript body .
therefore, all genes that exhibit length differences larger than the read length should be excluded from any of the reference gene sets .
reciprocal re-annotation reduces the number of false positive candidates
overview of the reciprocal re-annotation pipeline
to overcome problems due to length differences between orthologous genes and simultaneously maximize the number of comparable genes, we developed a pipeline to reciprocally re-annotate the reference genomes of the three species .
instead of directly annotating the d. simulans and d. mauritiana genomes individually using the d. melanogaster reference gene set, we first annotated the genome of d. simulans based on the protein set from d. melanogaster.
subsequently, we used these newly annotated d. simulans gene models to annotate the genome of d. mauritiana.
this gene set was then used as a reference to re-annotate again the previously generated d. simulans gene set.
and finally, we used these d. simulans gene models that already contain consensus features of d. simulans and d. mauritiana to re-annotate the d. melanogaster gene set .
therefore, we obtained the longest sequence present in all three species and then, if necessary, reduce its length in the other references accordingly.
thus, we expect to equalize the length of all the genes for the three references .
it is important to note here that it does not matter in which order the reciprocal re-annotation is done.
as long as the first reference is the one of d. melanogaster , we obtained the same results when we first annotate d. simulans or d. mauritiana .
reciprocal re-annotation efficiently reduces gene length differences between species
with the reciprocal re-annotation of the published genomes we obtained  <dig>  % of the  <dig>  d. melanogaster gene models in each of the three species .
in accordance with our expectations, only a small fraction of those genes found in all three species have a length difference of more than 49 bp :  <dig> genes  genes between d. melanogaster and d. mauritiana,  <dig> genes  between d. melanogaster and d. simulans and only  <dig> genes  between d. mauritiana and d. simulans.
hence, the reciprocal re-annotation of the published genomes allows the analysis of the highest number of comparable genes with less than 50 bp length differences in a differential gene expression study  of the  <dig>  d. melanogaster gene models).
evaluation of the reciprocal re-annotation in rna-seq experiments
to quantitatively test whether the number of false positives due to gene length differences is indeed reduced after reciprocal re-annotation, we applied a pairwise analysis of differential gene expression between d. melanogaster and d. mauritiana and d. simulans and d. mauritiana .
we mapped the rna-seq reads to the new references and performed the statistical analysis using the four methods evaluated above: deseq <dig>  deseq <dig> with length correction  <cit> , rpkm-limma  <cit> , rpkm-voom-limma  <cit> .
as for the published and directly re-annotated references, the statistical analysis with deseq <dig> resulted in the highest number of significantly differentially expressed genes .
this number clearly dropped to  <dig>  and  <dig>  % after including a length correction during the deseq <dig> analysis.
similarly, the number of differentially expressed genes is greatly reduced if rpkm-limma and rpkm-voom-limma are being used .
however, only  <dig>  and  <dig>  pairwise comparable genes exhibit length differences greater than 49 bp after reciprocal re-annotation.
one would expect that only those genes should be affected by any of the three length correction methods.
therefore, we propose that the combination of a reciprocal re-annotation in combination with a read-count based deseq <dig> analysis of differential gene expression is likely to provide the most comprehensive and reliable estimation of inter-specific gene expression differences.
this is further supported by the lack of a significant correlation between log2-fold changes and gene length differences if the deseq <dig> is used in combination with the reciprocal re-annotation as mapping reference .
although the correlation is not significant, we still find that most significantly differentially expressed genes with length differences larger than 49 bp have higher expression in the species with the longer transcript .
therefore, we propose that those genes should be filtered out from further differential gene expression analysis.
additionally, the seven genes that were validated using qpcr did not show a significant differential expression after their length was equalized , suggesting that the length correction during the annotation of genomes can facilitate the reduction in false positive candidate genes in rna-seq analyses.
assessment of power to detect differential gene expression
it is important to note that the gene models generated by our reciprocal re-annotation pipeline do not necessarily represent the complete gene and thus the most comprehensive annotation for each species.
this is due to the fact that potential full gene models in one species might have been adjusted to the shortest orthologous gene model.
therefore, in each round of annotation some gene models are truncated to fit the length of its orthologs .
if the gene models would be extremely shortened, this could of course lead to a loss of statistical power for the differential gene expression analysis.
in order to estimate how much sequence information we really lose, we compared the length of the d. melanogaster gene models before and after the reciprocal re-annotation.
this comparison shows that  <dig>   of the  <dig>  gene models still contain  <dig> to 100 % of their original sequence length after the reciprocal re-annotation .
next we assessed the potential loss of power by comparing the number of mapped reads between the published annotations  and the gene sets generated from our reciprocal re-annotation of the published genomes .
overall, the proportion of successfully mapped reads for all reference gene sets was between  <dig> and 67 % .
a large portion of this relatively low mapping rate can be explained by the fact that we excluded utr sequences from all reference gene sets, what accounts for about  <dig>  % of all mapped reads .
additionally, we only used the longest isoform of d. melanogaster for all annotations in the other two species .
therefore, some differentially spliced exons might not be represented in the newly generated gene sets.
however, the use of the sum of all exons only increases the mapping success by  <dig>  % if utrs are excluded and  <dig>  % if the utrs are included .
if the comparison of the expression of different isoforms across species is of interest one could perform the quantification on the level individual transcripts  <cit>  or even exons.
this approach requires of course a proper annotation of the different isoforms in all reference genomes and a dedicated mapping pipeline.
for our analysis, we found for all replicates more than  <dig> million mapped reads after reciprocal re-annotation  what has been shown to provide enough statistical power for differential gene expression analyses  <cit> .table  <dig> analysis of differential expression
log2fc
d. mel
d. mau
log2fc
log2fc
log2fc
log2fc
d. mel
d. mau
log2fc
expression comparison is for d. mauritiana vs. d. melanogaster, thus a positive log2-fold change  indicates higher expression in d. melanogaster and vice versa.
*p <  <dig> ; **p <  <dig> ; ***p <  <dig>
| list of rna-seq samples and the percentage and number of mapped reads to different reference transcriptomes
a reads originally 100 bp paired-end  were split in half to be 50 bp each and treated as single-end  reads
we observed an increase in the mapping percentage of up to 5 % in d. simulans and d. mauritiana when the reciprocally re-annotated gene sets are used as references .
this result shows that, although some gene models were now shorter, many genes that had been filtered out in the published genome annotations are actually expressed in these species.
the use of the re-annotated gene set only slightly decreases the mapping success by  <dig> to  <dig>  % in d. melanogaster , which is likely to be due to the artificial shortening of d. melanogaster gene models.
in summary, we show that the artificial shortening of transcripts after reciprocal re-annotation does not have a major impact on the power to detect differential gene expression.
practical considerations
we demonstrate that the use of all annotated exons instead of the longest isoform of each gene model does not significantly increase the power to detect differential gene expression.
in contrast, the inclusion of utr regions for the reciprocal re-annotation will clearly increase the number of mapped reads  and hence the statistical power.
however, the availability of utr sequence information strongly depends on the quality of the annotation of the species to compare, since utr and isoform predictions usually profit from the presence of rna-seq data to be incorporated in the annotation pipeline.
additionally, the annotation of utr regions might become more complicated if more distantly related species are studied, because utr regions tend to evolve faster than coding region  <cit> .
although we used very closely related species for our analysis here, we think that the presented reciprocal re-annotation is also applicable for genomes of more distantly related species.
as a consequence of a higher sequence divergence between distantly related species, inter-specific gene length differences are likely to be more pronounced.
if such genomes were used as mapping references, the direct use of length correction during the statistical analysis of differential gene expression might enhance the over-correction effect that we have demonstrated for the three presented methods.
additionally, if the gene lengths are very different between species, the length bias that has been reported for rpkm based normalization approaches  might be more pronounced.
therefore, we propose that the correction of the inter-specific length bias prior to read mapping using our reciprocal re-annotation pipeline should result in more robust results.
however, for more distantly related species, the reciprocal re-annotation is likely to result in more artificial shortening of the genes.
since this could reduce the power to detect differential gene expression, we propose to assess the length differences between species as we presented it here  prior to the sample preparation and sequencing and to adjust the coverage accordingly by generating more reads to increase sequencing depth.
in the presented case, at least one of the three drosophila species represents a well-established model system with a high quality genome assembly and annotation.
if this is available, the reciprocal re-annotation pipeline should of course be started with the highest quality annotation.
if the annotation quality of all genomes similar the pipeline could be started with any of the studied species, since we showed that the direction of the reciprocal annotation does not influence the final result.
however, if the quality of all annotations is comparably low, one should consider generating longer paired-end reads with higher coverage to first perform a de novo annotation with tools like augustus  <cit>  or braker <dig>  <cit>  using those reads to train the respective algorithm.
subsequently, the generated rna-seq reads can be used to assess differential gene expression using the reciprocally re-annotated references with length adjusted orthologous genes.
CONCLUSIONS
we have carried out a comprehensive comparison of the annotations of published genomes for the three closely related drosophila species, d. melanogaster, d. simulans and d. mauritiana.
this analysis reveals that different assembly strategies, annotation pipelines and filtering steps result in only a small fraction of genes that are comparable among all three species.
a direct re-annotation of the d. simulans and d. mauritiana genomes using the same d. melanogaster reference gene set and the same annotation pipeline significantly improves the comparability of the gene sets.
however, this direct re-annotation still results in length differences in many gene models between species.
based on these length differences between orthologous genes we tested four alternative methods to statistically assess differential gene expression using rna-seq, namely deseq <dig>  deseq <dig> with length correction, rpkm-limma and rpkm-voom-limma.
we show that none of these methods sufficiently accounts for the inter-specific gene length differences what is evident by a high number of false positive differentially expressed genes.
this finding is further supported by qpcr as an alternative transcript quantification method.
in order to further reduce the observed false positive rate, we argue that the length bias should be accounted for prior to the rna-seq analysis during the generation of the mapping references.
therefore, we implemented a robust reciprocal re-annotation pipeline that allows the generation of highly comparable gene sets to serve as mapping references for inter-specific rna-seq experiments.
applying rna-seq and qpcr we confirm the successful reduction of false positive candidate genes if the reciprocally re-annotated genomes are used as mapping references.
the reciprocal re-annotation pipeline can easily be adopted to re-annotate genomes of other closely related species or populations of animals and plants.
although we introduced our novel approach here to re-annotate three genomes at a time, it can of course be applied to two or more genomes.

 comparative studies of gene expression have been used to understand the regulation of a wide range of biological processes. with the development of next generation sequencing technologies, and in particular the use of illumina sequencing platforms, reliable genome wide comparison of gene expression between different biological conditions has become possible. moreover, a growing number of available genome and transcriptome sequences now provides the opportunity to compare gene expression not only in well-established, but also in emerging model systems. especially, the comparison of gene expression between both closely and distantly related species has great potential to help understand phenotypic divergence and species adaptations at a mechanistic level. experiments to study differential gene expression using ngs technologies are based on a sequencing library generated from reverse transcribed messenger rna that is extracted from the tissue and conditions of interest. illumina sequencing, for example, results in the generation of millions of short reads ranging from 36 bp to 150 bp. the first step of the bioinformatics analysis is to align these reads to a reference that represents all transcripts that should be quantified. this reference can be a whole genome sequence with annotated gene models or a transcriptome. the latter can either be generated by a de novo assembly of the rna-seq reads or it could be extracted from an annotated genome. the next step is to determine the number of reads that are aligned to a gene model or transcript. depending on the type of reference used various different methods have been established. finally, the number of reads assigned to a given gene model or transcript is compared between different conditions to identify differentially expressed genes. the steps outlined above for a general rna-seq experiment are suitable to compare gene expression levels between different conditions, stages or tissues of the same species. however, comparison of gene expression between different species or populations of the same species needs to account for differences in gene sequences. in this case, reads should be mapped to species-specific references for which the expression level of a gene in one of the species is compared to the expression level of its ortholog in the other species. most importantly, this requires sets of orthologous genes reliably identified in all references. since genomes or transcriptomes are usually generated by different research groups for different applications and using different pipelines for assembly and annotation, annotated references for inter-specific gene expression studies are often not comparable. for instance, orthologous genes might be missing from one or more of the references as result of natural variation or technical problems like incomplete assemblies or too many sequencing errors, which hampers unequivocal identification of orthologous genes. additionally, it is common practice to filter out genes that are incomplete or lack synteny in relation to a model reference from new gene model predictions. even though there are many tools available to perform genome annotation, a general standard does not exist. therefore, the final gene set generated by each genome project will have genes missing as a result of methodological problems and filtering criteria, and this can directly influence the result of the differential gene expression analysis. additionally, even if most one-to-one orthologs have been successfully identified in different references, these gene models may vary in length for various reasons: first, the genes could naturally differ in length among species. second, as a consequence of differences in the sequence or assembly quality of the reference genomes, orthologous gene models might be truncated in one or more of the references. to our knowledge, a comprehensive evaluation of methods that could be applied to account for inter-specific gene length differences has not been performed yet. a plethora of statistical approaches have been developed to determine whether differences in the number of reads are due to technical variation or due to real biological differences in gene expression. detailed evaluation and comparison of these methods concluded that the most accurate statistical validation of differential gene expression is reached when statistical models are used that directly take the number of aligned reads into account. these methods include standard and generalized poisson and negative binomial distributions to model count-based expression data as implemented in deseq, deseq, edger or degps. also the differential expression analysis based on moderated t-statistics as implemented in the limma package using log-transformed count per million values originating from normalization with voom has been shown to perform extremely well. while all of these methods account for most technical biases and control well false positive rates, none of these methods is specifically designed to account for gene length differences as they occur in inter-specific expression studies. one potential solution could be the application of the normalization method reads per kilobase per million mapped reads as it accounts for length differences in gene models. however, it has been shown that even after correcting for length differences, a longer transcript is more likely to appear as differentially expressed if rpkm values are used to assess the statistical significance. rpkm normalization is still widely used to compare gene expression levels of different genes within a species, but to our knowledge it has not been tested if this method efficiently normalizes length differences when comparing gene expression in different species. here we show that the published genomes of three closely related drosophila species, d. melanogaster, d. simulans and d. mauritiana have qualitative limitations as references for comparative gene expression studies. this is mainly due to the fact that many genes cannot be properly compared because orthologous genes are missing in the annotation of at least one of the genomes. even after a direct re-annotation of the three genomes using the same annotation pipeline many orthologous gene models exhibit significant length differences. taking advantage of these inter-specific gene length differences in the published and the directly re-annotated references, we benchmarked four statistical frameworks for their ability to reduce the number of potentially false positives. we demonstrate that none of these methods sufficiently accounts for the observed differences in gene length. therefore, we propose that the length normalization should be performed prior to read mapping during the generation of the mapping references. we report a straightforward re-annotation method that relies on a reciprocal re-annotation of orthologous gene models in two or more species. this approach allows the comparison of nearly all genes that have been annotated in d. melanogaster in all three species. additionally, we find that the use of the new gene sets as mapping references results in a more robust estimation of transcript abundance and a more reliable comparison of gene expression levels between species. we propose that the generation and annotation of new genome resources or the re-annotation of existing genomes will be powerful tools to establish gene expression profiling in many emerging model systems.  analysis of published genome annotations reveals a reduced number of comparable gene models for differential gene expression studies between species we first assessed the completeness and comparability of the published gene sets for the three closely related species d. melanogaster, d. simulans and d. mauritiana. at the time of our analysis, the annotation of the d. melanogaster genome - one of the best curated metazoan genomes available at flybase - included unique protein coding genes. the most recent annotations for d. simulans and d. mauritiana were generated using the d. melanogaster gene set as reference. both gene sets contain a large fraction of the d. melanogaster genes. however, orthologs of almost d. melanogaster genes are not included in each of the final gene sets either because the authors applied various filtering steps to exclude incomplete orthologous sequence with respect to the d. melanogaster gene or because the genes are not present in one of the species. since these filtering steps are influenced by the quality of each of the assembled genome and the scientific question of each research group, the missing genes in both annotations are not the same. only genes can be identified unequivocally as orthologs in all three annotations. among the genes missing in at least one annotation, we found some important and well-studied developmental genes including the hox genes abdominal b, ultrabithorax or antennapedia, the head and brain patterning gene orthodenticle and the segment polarity gene hedgehog.table number of genes obtained by each annotation method d. melanogaster d. simulans d. mauritiana the last column contains the number of genes for which 1: orthologs were identified in the three species. “after filtering” indicates the remaining common genes after filtering out genes with length difference larger than bp. percentages in brackets are always given in relation to the total number of gene models in d. melanogaster next we assessed the comparability of the three reference genome annotations in terms of gene length, since length differences larger than the length of the rna-seq reads are likely to introduce a bias during mapping and the subsequent differential expression analysis. if we consider 50 bp single-end reads, which have been shown to be long enough to produce accurate results when measuring differential gene expression, genes that have a length difference larger than 49 bp among the annotations of the three drosophila species are likely to bias a subsequent differential gene expression analysis. a pairwise comparison of annotated gene length for the genes present in all three drosophila species shows that in the published annotations, the gene length differences are larger than 49 bp in  % of the orthologous genes between d. mauritiana and d. simulans,  % between d. melanogaster and d. simulans and  % between d. melanogaster and d. mauritiana. if these length differences are not accounted for, these genes could result in false positives in a differential gene expression analysis.fig. pairwise length difference between orthologous genes. bars indicate the number of genes with that difference in length for each pair of species. green shades indicate differences lower than 50 bp while orange to red indicate larger differences. the comparison is showed for the published annotations, b the direct re-annotation of the published genomes and c the reciprocal re-annotation of the published genomes direct re-annotation of published genomes next we asked whether a direct re-annotation of the d. simulans and d. mauritiana genomes individually using the same d. melanogaster gene set as reference and the same annotation pipeline allows the comparison of more genes in an inter-specific differential gene expression study. we used the protein sequences of d. melanogaster as reference to re-annotate the published genomes of d. simulans and d. mauritiana using the program exonerate. without applying any filtering, we find orthologs of d. melanogaster genes that are comparable among the two species. next, we determined the length of directly re-annotated genes that are found in all three species. this comparison shows an increase in number of genes with the same length between the three species after direct re-annotation. however, a high number of orthologous genes have a length difference of more than 49 bp: between d. melanogaster and d. mauritiana, between d. melanogaster and d. simulans and between d. mauritiana and d. simulans. these observed length differences could be due to real natural variation in coding sequence length between species or they could be technical artifacts, for example truncated gene sequences arising from sequencing or genome assembly errors. in summary, although annotated genomes are available for the three closely related drosophila species d. melanogaster, d. simulans and d. mauritiana, their use as mapping references for inter-species differential gene expression analyses is limited due to missing orthologs and a potential bias because of different annotated gene lengths. the use of the same reference gene set, annotation pipeline and the lack of filtering incomplete gene sequences results in an increase in the number of comparable genes in these three closely related species. length difference in reference genes introduces biases in differential expression studies since we find a high number of gene models with length differences > 49 bp in the published annotations and after the direct re-annotation, the three drosophila genomes are excellent models to test whether length differences larger than the read length do indeed influence the statistical analysis of differential gene expression. we used the published d. melanogaster annotation and the newly generated direct re-annotations of d. simulans and d. mauritiana as mapping references for pairwise comparisons of gene expression between d. melanogaster and d. mauritiana and d. simulans and d. mauritiana using 50 bp single-end illumina rna-seq reads generated for these three species. rna-seq reads generated from one species were mapped only to the gene set of that species. using this experimental setup, we compared four different statistical frameworks, namely deseq deseq with gene length correction, limma with length correction based on rpkm and voom-limma including rpkm length correction. for each method, we first report the number of differentially expressed genes for each of the two pairwise species comparisons. next we evaluate the impact of the length differences between gene models on the fold-change in gene expression between species. and eventually, we compare the results of each of the four methods to an independent qpcr experiment for a subset of genes. deseq without length correction first we performed the statistical analysis for differential gene expression with deseq using directly the read counts for each gene model. for both pairwise comparisons using the published and the direct re-annotation as reference, we found that between and  % of all comparable genes are significantly differentially expressed. additionally, we found a very strong correlation between inter-specific length differences of the gene models and the log2-fold change in gene expression. the negative correlation means that genes that are longer in one species appear to be more up-regulated and vice versa. the correlation can be explained by the mapping procedure: in orthologous genes with different length, more reads align to the ortholog that has the longer gene model. this results in an artificially higher expression for this specific gene in the species with the longer gene model. from this correlation we also see that most of those genes with length differences and a high log2-fold change are also significantly differentially expressed, showing that this method introduces a large number of false positives.fig. length differences between orthologous genes introduce gene expression biases. relation between length differences and the log2-fold change in the rna-seq experiment between d. mauritiana and d. simulans using the direct re-annotation of these species as mapping references. dots represent genes with length difference > 49 bp in these annotations. genes significantly differentially expressed in the presented analysis are shown in red. a negative log2-fold change indicates higher expression in d. mauritiana. a positive length difference indicates that the ortholog of d. mauritiana is longer. the p-value and rho of the spearman’s rank correlation are indicated on the lower right side of the plots. a results of deseq without length correction. b results of deseq applying length normalization factor matrix. c results of applying rpkm normalization and limma to call differentially expressed genes. d results of applying voom normalization followed by a length normalization matrix and limma to call differentially expressed genes fig. schematic representation of length bias in inter-species differential expression analysis and our reciprocal re-annotation strategy to correct it. a length bias in the analysis of a non-differentially expressed gene. coloured rectangles represent the part of the transcript which is included as reference for the rna-seq reads to map to, while unfilled rectangles are regions of the transcript which are omitted and to which rna-seq reads cannot be mapped. red “n”s represent sequencing errors that prevent the complete annotation of a transcript. mapped reads are shown as thin black lines and the number bellow indicates the total of reads mapped. if one transcript is shorter in one of the references compared to its orthologs, for the same expression levels fewer reads will map to it. this can result in false positives in the analysis of differential expression. our strategy to correct this bias is to shorten the orthologs in the other references to match the length of the shorter sequence. b pipeline of reciprocal transcriptome re-annotation method. black numbers in white circles represent genome annotation steps using the “est2genome” command of exonerate. grey numbers in grey circles represent conversion of the resulting gff file into a new transcript set. filled horizontal bars represent the annotated set of transcripts; non-filled horizontal bars at the start/end of the transcripts represent parts of the transcript that cannot be correctly annotated in one reference and are therefore eliminated from the transcript set. the boxes with red frame indicate the transcript sets that will be used as reference for rna-seq read mapping. step 1: the transcript set of the best annotated genomes is used to annotate one of the other genomes and generate a new transcript set for this species. due to sequencing errors, some transcripts will be shorter. step 2: the new transcript set form d. simulans is used to annotate the last genome. the gene set generated contains shorter transcripts due to sequencing errors in d. mauritiana but also in d. simulans. step 3: the transcript set from d. mauritiana is used to re-annotate the previously generated set from d. simulans to integrate the information from the d. mauritiana assembly. step 4: the second transcript set from d. simulans is used to annotate the d. melanogaster set in order to integrate the information from d. simulans and d. mauritiana in order to specifically test whether differences in the length of gene models indeed influence the differential expression analysis we chose seven genes that were shorter in the d. mauritiana published annotation compared to the published d. melanogaster annotation. when we analysed the differential expression using deseq without any length correction, the expression of all seven genes were significantly different. the log2-fold change value indicated that d. melanogaster had a significantly higher expression than d. mauritiana. to validate the results obtained by the rna-seq experiment, we measured the relative expression of the seven genes in d. melanogaster and d. mauritiana using qpcr. this analysis showed that the seven genes that had length differences in the species-specific annotations were not significantly differentially expressed. as a control we chose another three genes that showed significant differential expression in the rna-seq data but had the same length in both species in the two annotation methods. we found that piwi and alrm showed a significantly higher expression in d. mauritiana when using this alternative quantification method, confirming the results obtained by rna-seq. moreover, nplp had higher expression in d. melanogaster again consistent with our rna-seq data, although this difference was not significant.table differentially expressed genes and correlation between calculated log2-fold changes and length difference between orthologous genes. results are shown for the four applied methods, the three studied annotation strategies and the two described pairwise species comparisons d. mau vs d. mel d. mau vs d. sim d. mau vs d. mel d. mau vs d. sim d. mau vs d. mel d. mau vs d. sim d. mau vs d. mel d. mau vs d. sim a d. mau vs d. mel a d. mau vs d. sim d. mau vs d. mel d. mau vs d. sim d. mau vs d. mel d. mau vs d. sim d. mau vs d. mel a d. mau vs d. sim d. mau vs d. mel d. mau vs d. sim d. mau vs d. mel d. mau vs d. sim a d. mau vs d. mel a d. mau vs d. sim a d. mau vs d. mel d. mau vs d. sim afdr bspearman’s rank correlation is measured between log2fc and length difference of genes showing more than 49 bp length difference published annotation: d. mau vs. d.mel: genes/d.mau vs. d.sim: genes; direct annotation: d.mau vs. d.mel: genes/d.mau vs. d.sim: genes; reciprocal annotation: d.mau vs. d.mel: genes/d.mau vs. d.sim: genes *p < ; **p < ; ***p <  fig. qpcr results. boxplot of normalized ct values for each studied gene boxplot is showed for ct values in d. melanogaster orer and d. mauritiana tam. or t-welch-test; *p < , **p < ; ***p <  ) in summary, these results suggest a high level of potentially false positive candidates when methods based on direct read counts without the application of length correction are used with mapping references that exhibit differences in the length of orthologous genes. deseq with length correction next we benchmarked the use of deseq including a normalization factor matrix incorporating gene length to account for the length differences between orthologous genes. using this approach for pairwise gene expression analyses, we found that only to  % of the comparable genes were significantly differentially expressed. even though the gene length was accounted for during the deseq analysis of differential gene expression, we still find a correlation between inter-specific gene length differences and log2-fold changes. however, the significance of this correlation is greatly reduced in comparison to the deseq analysis without length correction, suggesting that the length correction incorporated in deseq helps to reduce the number of false positive candidate genes. this finding is further supported by the comparison of the rna-seq results to the qpcr data. after length correction only one of the seven genes that are longer in d. melanogaster show significant differential expression. limma with rpkm length correction rpkm values are commonly calculated for rna-seq datasets to account for variation in library sizes and to correct for length differences between different genes within the same species. the moderated t-statistics incorporated in the limma r package can subsequently be used to assess differential gene expression. it has not been tested if this approach also corrects properly for differences in the length of the same gene being compared between two species. using this method, we found between  % and to  % of the comparable genes to be significantly differentially expressed. our correlation analysis shows that the correction of a length bias using rpkm values still results in a clearly significant correlation between the gene length difference and the observed log2-fold change. however, compared to the deseq analysis without length correction, the significance values are highly reduced, showing that the number of false positives is lower. accordingly, six of the seven genes that we benchmarked with qpcr show no significant differential gene expression although they show clear length differences between d. melanogaster and d. mauritiana. again cp is the only gene that appears as significantly differentially expressed also after correcting for length differences. voom-limma with rpkm length correction it has recently been shown that differential gene expression analysis with limma using normalized read counts from voom perform very well for rna-seq datasets. although this method is designed to work with direct read counts, in this case we tested it with an additional transcript length correction. between and  % of the comparable genes are significantly differentially expressed. after length correction and normalization with voom, we found a significant correlation between gene length differences and log2-fold changes when the published annotations and the directly re-annotated reference gene sets were used however, this was slightly reduced compared to the rpkm-limma analysis, especially for the d. simulans and d. mauritiana comparison.. for the seven qpcr benchmarked genes that have clear length differences between d. melanogaster and d. mauritiana the use of the voom-limma method results in three significantly differentially expressed genes, suggesting a higher false positive rate. length correction during the statistical analysis might be insufficient the comprehensive comparison of four methods for differential gene expression analysis shows that the incorporation of a length correction drastically reduces the number of false positive candidate genes. although the correlation between gene length differences and the observed log2-fold changes is reduced in the three methods that account for gene length differences, none of them sufficiently corrects the length bias present in the two gene sets used as mapping references. this is also supported by the qpcr validation of seven genes that exhibit clear length differences between the published d. melanogaster and d. mauritiana annotations. in all three methods at least one gene was still artificially significantly differentially expressed. this is most pronounced for the voom-limma method where three of the seven genes are significantly differentially expressed. of the seven genes we analysed using qpcr, cp was in all cases identified as a false positive candidate. in order to further characterise this gene, we visually inspected the distribution of mapped reads. interestingly, the 3′-region is missing in the d. mauritiana ortholog of cp and we found clearly more d. melanogaster reads that map to this 3′-part of the transcript than to the 5′-region. the number of d. melanogaster and d. mauritiana reads mapped to the 5′-part of the transcript is comparable. hence, a very likely explanation for the inefficient length correction of the three applied methods could be an unequal distribution of the mapped reads along the transcripts. besides an insufficient length correction, the deseq method including a length matrix results in the lowest number of significantly differentially expressed genes, suggesting that the length correction applied here might be extremely conservative and could lead to a high rate of false negatives. interestingly, in many pairwise comparisons we found more significantly differentially expressed genes with a higher expression in d. mauritiana compared to d. melanogaster and d. simulans, although d. mauritiana gene models are generally shorter than those of the other two species. this finding suggests that the length correction applied here might reduce the power to detect differential expression for the already short d. mauritiana genes. in summary, all three methods that include a length correction decrease the chance of identifying false positives. the rpkm-voom-limma and rpkm-limma methods seem to give the best ratio of false positives and false negatives, while deseq including a length matrix is very conservative. however, none of the length correction methods tested does efficiently account for all differences in gene length observed in the reference annotation of the three studied drosophila species. the length bias is most obvious when the distribution of reads is not uniform across the transcript body. therefore, all genes that exhibit length differences larger than the read length should be excluded from any of the reference gene sets. reciprocal re-annotation reduces the number of false positive candidates overview of the reciprocal re-annotation pipeline to overcome problems due to length differences between orthologous genes and simultaneously maximize the number of comparable genes, we developed a pipeline to reciprocally re-annotate the reference genomes of the three species. instead of directly annotating the d. simulans and d. mauritiana genomes individually using the d. melanogaster reference gene set, we first annotated the genome of d. simulans based on the protein set from d. melanogaster. subsequently, we used these newly annotated d. simulans gene models to annotate the genome of d. mauritiana. this gene set was then used as a reference to re-annotate again the previously generated d. simulans gene set. and finally, we used these d. simulans gene models that already contain consensus features of d. simulans and d. mauritiana to re-annotate the d. melanogaster gene set. therefore, we obtained the longest sequence present in all three species and then, if necessary, reduce its length in the other references accordingly. thus, we expect to equalize the length of all the genes for the three references. it is important to note here that it does not matter in which order the reciprocal re-annotation is done. as long as the first reference is the one of d. melanogaster, we obtained the same results when we first annotate d. simulans or d. mauritiana. reciprocal re-annotation efficiently reduces gene length differences between species with the reciprocal re-annotation of the published genomes we obtained  % of the d. melanogaster gene models in each of the three species. in accordance with our expectations, only a small fraction of those genes found in all three species have a length difference of more than 49 bp: genes genes between d. melanogaster and d. mauritiana, genes between d. melanogaster and d. simulans and only genes between d. mauritiana and d. simulans. hence, the reciprocal re-annotation of the published genomes allows the analysis of the highest number of comparable genes with less than 50 bp length differences in a differential gene expression study of the d. melanogaster gene models). evaluation of the reciprocal re-annotation in rna-seq experiments to quantitatively test whether the number of false positives due to gene length differences is indeed reduced after reciprocal re-annotation, we applied a pairwise analysis of differential gene expression between d. melanogaster and d. mauritiana and d. simulans and d. mauritiana. we mapped the rna-seq reads to the new references and performed the statistical analysis using the four methods evaluated above: deseq deseq with length correction, rpkm-limma, rpkm-voom-limma. as for the published and directly re-annotated references, the statistical analysis with deseq resulted in the highest number of significantly differentially expressed genes. this number clearly dropped to and  % after including a length correction during the deseq analysis. similarly, the number of differentially expressed genes is greatly reduced if rpkm-limma and rpkm-voom-limma are being used. however, only and pairwise comparable genes exhibit length differences greater than 49 bp after reciprocal re-annotation. one would expect that only those genes should be affected by any of the three length correction methods. therefore, we propose that the combination of a reciprocal re-annotation in combination with a read-count based deseq analysis of differential gene expression is likely to provide the most comprehensive and reliable estimation of inter-specific gene expression differences. this is further supported by the lack of a significant correlation between log2-fold changes and gene length differences if the deseq is used in combination with the reciprocal re-annotation as mapping reference. although the correlation is not significant, we still find that most significantly differentially expressed genes with length differences larger than 49 bp have higher expression in the species with the longer transcript. therefore, we propose that those genes should be filtered out from further differential gene expression analysis. additionally, the seven genes that were validated using qpcr did not show a significant differential expression after their length was equalized, suggesting that the length correction during the annotation of genomes can facilitate the reduction in false positive candidate genes in rna-seq analyses. assessment of power to detect differential gene expression it is important to note that the gene models generated by our reciprocal re-annotation pipeline do not necessarily represent the complete gene and thus the most comprehensive annotation for each species. this is due to the fact that potential full gene models in one species might have been adjusted to the shortest orthologous gene model. therefore, in each round of annotation some gene models are truncated to fit the length of its orthologs. if the gene models would be extremely shortened, this could of course lead to a loss of statistical power for the differential gene expression analysis. in order to estimate how much sequence information we really lose, we compared the length of the d. melanogaster gene models before and after the reciprocal re-annotation. this comparison shows that of the gene models still contain to 100 % of their original sequence length after the reciprocal re-annotation. next we assessed the potential loss of power by comparing the number of mapped reads between the published annotations and the gene sets generated from our reciprocal re-annotation of the published genomes. overall, the proportion of successfully mapped reads for all reference gene sets was between and 67 %. a large portion of this relatively low mapping rate can be explained by the fact that we excluded utr sequences from all reference gene sets, what accounts for about  % of all mapped reads. additionally, we only used the longest isoform of d. melanogaster for all annotations in the other two species. therefore, some differentially spliced exons might not be represented in the newly generated gene sets. however, the use of the sum of all exons only increases the mapping success by  % if utrs are excluded and  % if the utrs are included. if the comparison of the expression of different isoforms across species is of interest one could perform the quantification on the level individual transcripts or even exons. this approach requires of course a proper annotation of the different isoforms in all reference genomes and a dedicated mapping pipeline. for our analysis, we found for all replicates more than million mapped reads after reciprocal re-annotation what has been shown to provide enough statistical power for differential gene expression analyses.table analysis of differential expression log2fc d. mel d. mau log2fc log2fc log2fc log2fc d. mel d. mau log2fc expression comparison is for d. mauritiana vs. d. melanogaster, thus a positive log2-fold change indicates higher expression in d. melanogaster and vice versa. *p < ; **p < ; ***p <  | list of rna-seq samples and the percentage and number of mapped reads to different reference transcriptomes a reads originally 100 bp paired-end were split in half to be 50 bp each and treated as single-end reads we observed an increase in the mapping percentage of up to 5 % in d. simulans and d. mauritiana when the reciprocally re-annotated gene sets are used as references. this result shows that, although some gene models were now shorter, many genes that had been filtered out in the published genome annotations are actually expressed in these species. the use of the re-annotated gene set only slightly decreases the mapping success by to  % in d. melanogaster, which is likely to be due to the artificial shortening of d. melanogaster gene models. in summary, we show that the artificial shortening of transcripts after reciprocal re-annotation does not have a major impact on the power to detect differential gene expression. practical considerations we demonstrate that the use of all annotated exons instead of the longest isoform of each gene model does not significantly increase the power to detect differential gene expression. in contrast, the inclusion of utr regions for the reciprocal re-annotation will clearly increase the number of mapped reads and hence the statistical power. however, the availability of utr sequence information strongly depends on the quality of the annotation of the species to compare, since utr and isoform predictions usually profit from the presence of rna-seq data to be incorporated in the annotation pipeline. additionally, the annotation of utr regions might become more complicated if more distantly related species are studied, because utr regions tend to evolve faster than coding region. although we used very closely related species for our analysis here, we think that the presented reciprocal re-annotation is also applicable for genomes of more distantly related species. as a consequence of a higher sequence divergence between distantly related species, inter-specific gene length differences are likely to be more pronounced. if such genomes were used as mapping references, the direct use of length correction during the statistical analysis of differential gene expression might enhance the over-correction effect that we have demonstrated for the three presented methods. additionally, if the gene lengths are very different between species, the length bias that has been reported for rpkm based normalization approaches might be more pronounced. therefore, we propose that the correction of the inter-specific length bias prior to read mapping using our reciprocal re-annotation pipeline should result in more robust results. however, for more distantly related species, the reciprocal re-annotation is likely to result in more artificial shortening of the genes. since this could reduce the power to detect differential gene expression, we propose to assess the length differences between species as we presented it here prior to the sample preparation and sequencing and to adjust the coverage accordingly by generating more reads to increase sequencing depth. in the presented case, at least one of the three drosophila species represents a well-established model system with a high quality genome assembly and annotation. if this is available, the reciprocal re-annotation pipeline should of course be started with the highest quality annotation. if the annotation quality of all genomes similar the pipeline could be started with any of the studied species, since we showed that the direction of the reciprocal annotation does not influence the final result. however, if the quality of all annotations is comparably low, one should consider generating longer paired-end reads with higher coverage to first perform a de novo annotation with tools like augustus or braker using those reads to train the respective algorithm. subsequently, the generated rna-seq reads can be used to assess differential gene expression using the reciprocally re-annotated references with length adjusted orthologous genes.  we have carried out a comprehensive comparison of the annotations of published genomes for the three closely related drosophila species, d. melanogaster, d. simulans and d. mauritiana. this analysis reveals that different assembly strategies, annotation pipelines and filtering steps result in only a small fraction of genes that are comparable among all three species. a direct re-annotation of the d. simulans and d. mauritiana genomes using the same d. melanogaster reference gene set and the same annotation pipeline significantly improves the comparability of the gene sets. however, this direct re-annotation still results in length differences in many gene models between species. based on these length differences between orthologous genes we tested four alternative methods to statistically assess differential gene expression using rna-seq, namely deseq deseq with length correction, rpkm-limma and rpkm-voom-limma. we show that none of these methods sufficiently accounts for the inter-specific gene length differences what is evident by a high number of false positive differentially expressed genes. this finding is further supported by qpcr as an alternative transcript quantification method. in order to further reduce the observed false positive rate, we argue that the length bias should be accounted for prior to the rna-seq analysis during the generation of the mapping references. therefore, we implemented a robust reciprocal re-annotation pipeline that allows the generation of highly comparable gene sets to serve as mapping references for inter-specific rna-seq experiments. applying rna-seq and qpcr we confirm the successful reduction of false positive candidate genes if the reciprocally re-annotated genomes are used as mapping references. the reciprocal re-annotation pipeline can easily be adopted to re-annotate genomes of other closely related species or populations of animals and plants. although we introduced our novel approach here to re-annotate three genomes at a time, it can of course be applied to two or more genomes. 
a comprehensive evaluation of four statistical frameworks shows that none of these methods sufficiently accounts for inter-specific gene length differences, which inevitably results in false positive candidate genes. we show that the published genomes and annotations of the three closely related drosophila species d. melanogaster, d. simulans and d. mauritiana have limitations for inter-specific gene expression studies. this is due to missing gene models in at least one of the genome annotations, unclear orthology assignments and significant gene length differences in the different species. we present a straight-forward reciprocal re-annotation pipeline that allows to reliably compare the expression for nearly all genes annotated in d. melanogaster. we conclude that our reciprocal re-annotation of previously published genomes facilitates the analysis of significantly more genes in an inter-specific differential gene expression study. until now, the main applications of this technique have been to study the variation of gene expression in a whole organism, tissue or cell type under different conditions or at different developmental stages. we propose that published reference genomes should be re-annotated before using them as references for rna-seq experiments to include as many genes as possible and to account for a potential length bias. rna-seq based on short reads generated by next generation sequencing technologies has become the main approach to study differential gene expression. however, rna-seq also has a great potential to be used in evolutionary studies to investigate gene expression divergence in closely related species. rna-seqannotationdifferential gene expressionexoneratedrosophilaclosely related speciesemerging model systemsrpkmdeseq2voomlimmalength biashttp://dx.doi.org/ /501100001663volkswagen foundation 983posnien nico http://dx.doi.org/ /501100000781european research council242553mcgregor alistair p. issue-copyright-statement© the author 2016
 rna-seq based on short reads generated by next generation sequencing technologies has become the main approach to study differential gene expression. until now, the main applications of this technique have been to study the variation of gene expression in a whole organism, tissue or cell type under different conditions or at different developmental stages. however, rna-seq also has a great potential to be used in evolutionary studies to investigate gene expression divergence in closely related species.  we show that the published genomes and annotations of the three closely related drosophila species d. melanogaster, d. simulans and d. mauritiana have limitations for inter-specific gene expression studies. this is due to missing gene models in at least one of the genome annotations, unclear orthology assignments and significant gene length differences in the different species. a comprehensive evaluation of four statistical frameworks shows that none of these methods sufficiently accounts for inter-specific gene length differences, which inevitably results in false positive candidate genes. we propose that published reference genomes should be re-annotated before using them as references for rna-seq experiments to include as many genes as possible and to account for a potential length bias. we present a straight-forward reciprocal re-annotation pipeline that allows to reliably compare the expression for nearly all genes annotated in d. melanogaster.  we conclude that our reciprocal re-annotation of previously published genomes facilitates the analysis of significantly more genes in an inter-specific differential gene expression study. we propose that the established pipeline can easily be applied to re-annotate other genomes of closely related animals and plants to improve comparative expression analyses. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords rna-seqannotationdifferential gene expressionexoneratedrosophilaclosely related speciesemerging model systemsrpkmdeseq2voomlimmalength biashttp://dx.doi.org/ /501100001663volkswagen foundation 983posnien nico http://dx.doi.org/ /501100000781european research council242553mcgregor alistair p. issue-copyright-statement© the author 2016 