BACKGROUND
rna-seq  has the potential to transform the way we study gene structure and expression  <cit> .
the ability to identify novel exons and splice sites  <cit>  is just the beginning.
although there have been claims that rna-seq is more sensitive  <cit>  and has a larger dynamic range  <cit>  than a microarray, these claims are now being challenged  <cit> .
for all the optimism surrounding rna-seq, there is growing evidence that estimating gene expression in an rna-seq environment is not straightforward.
in addition, technical variance, sequencing bias and mapping bias have been reported  <cit> .
the power and promise of rna-seq technology demand our attention.
new papers on experimental design  <cit>  will lead the way to more thoughtful experimentation.
one important component in careful study planning is an understanding of the different sources of variability so they can be accounted for in the experimental design  <cit> .
variability can occur at many levels, both biological and technical.
biological variability is unaffected by the technology, as genotypes, individuals and even cell types vary regardless of how they are measured.
bullard et al.
<cit>  examined the effects of library construction, flow cell and lane on detection of differential expression.
the data used there was commercial grade rna and the samples compared expected to be very divergent.
they concluded that the biological variation was larger than technical variation, underscoring the importance of including biological replicates in the study design.
the observation that biological variability is large does not imply that technical variability should be ignored in the experimental design of a study.
high levels of noise in a technology must be considered during experimental design  <cit> .
marioni et al.
<cit>  demonstrate clearly that technical variability is not different from what is expected due to random noise .
yet, the magnitude of the sampling error is not explicitly discussed.
one particular feature of deep sequencing is that even though there are millions of reads, the proportion of the mrna that is actually sequenced is low.
given a solexa/illumina rna-seq library of  <dig> nm concentration with a mean insert size of  <dig> bp and a volume of  <dig> ul, the number of molecules in a typical rna library is estimated to be  <dig>  ×  <dig>  on the solexa/illumina gaiix, approximately  <dig> million of the total possible molecules are sampled in a given lane.
this represents approximately  <dig> % of the total number of available molecules.
even with the anticipated increase in the next generation of technology, the sampling fraction is will be less than  <dig> %.
the impact of such a low sampling fraction on the statistical properties of rna-seq is the focus of this manuscript.
we address the following questions: 1) does a substantial amount of technical variability exist?
2) is the impact of technical variability the same for all levels of coverage?
the study presented here considers technical variation in three experiments using solexa/illumina technology.
although the number of exons consistently detected is improved by increasing the number of mappable reads, the detection of exons with low levels of coverage is inconsistent among technical replicates.
inconsistent detection of exons is most pronounced for exons with average coverage of less than  <dig> reads per nucleotide.
in addition, there can be substantial disagreement in the estimated level of expression among technical replicates, even when the coverage is high.
although the technical variability is not unexpected, technical variability cannot be ignored in rna-seq experiments.
methods
experiments
experiment 1: three independent samples of d. melanogaster female heads were collected with each sample representing a unique pool of biological material.
each sample was prepared according to manufacturer's instructions and then the same library was run on two lanes of a solexa/illumina flow cell, resulting in two technical replicates for each biological replicate, runs were  <dig> base-pair paired end.
experiment 2: three independent samples of d. simulans male heads were collected with each sample representing a unique pool of biological material.
each sample was prepared according to manufacturer's instructions and then the same library was run on two lanes of a solexa/illumina flow cell, resulting in two technical replicates for each biological replicate, runs were  <dig> base-pair paired end.
experiment 3: one sample of d. melanogaster cell line c <dig> was run on  <dig> lanes yielding  <dig> technical replicates for  <dig> biological replicate   <cit> .
data from these experiments were  <dig> bp reads.
the experimental design for experiments  <dig> -  <dig> are outlined in figure  <dig>  all samples were run using solexa/illumina paired end procedures.
the relationship of lanes  is unknown for all three experiments.
the first two experiments were conducted by joe dunham and michele arbeitman at usc and sequenced at oregon state core facility while the third experiment was conducted by susan celniker at lawrence livermore national laboratory .
simulation
the number of molecules present in the library is estimated to be  <dig>  ×  <dig>  this estimate was derived using the illumina protocol.
specifically, the starting material was assumed to be  <dig> ng of mrna that then resulted in  <dig> ng of library with an insert size of 250bp in a volume of  <dig> ul for a final library concentration of  <dig> nm.
the number of pmoles is therefore  <dig> .
the number of molecules is found by multiplying avogadro's number which is molecules/moles and adjusting for units  to give  <dig>  ×  <dig> molecules.
on the current solexa/illumina technology, the gaiix, approximately  <dig> million of the total possible molecules are sampled in a given lane.
this represents approximately  <dig> %  of the total number of available molecules for analysis .
to understand how random sampling may produce variation in technical replicates in an rna-seq experiment, a small simulation of two technical replicates was conducted.
we assume that the approximately  <dig>  ×  <dig> molecules of mrna were transcribed from  <dig>  genes.
gene expression for the  <dig>  genes was modeled using a gamma density function with shape parameter  <dig> and scale parameter  <dig>  this ensures that many genes will have a small number of molecules assigned to them and some genes will have a large number of molecules assigned to them.
the range of the gamma density function was examined over the interval from  <dig> to  <dig>  the integers from  <dig> to  <dig> represent bins of genes with similar relative levels of expression.
to simulate relative expression for each of the i =  <dig>   <dig>  ...,  <dig>  individual genes, a random integer  from a discrete uniform distribution between  <dig> and  <dig> were drawn and the value of the gamma density ti was calculated at si.
the number of total molecules of mrna assigned to gene i, gi, was, ni= ti/Σ *  <dig>  ×  <dig> rounded to the nearest integer.
Σ = z. due to the rounding, z was approximately, but not exactly equal to <dig>  ×  <dig>  for each gi with ni molecules a unique integer between  <dig> and z was assigned, resulting in a unique association between an integer and a gene.
to simulate the sampling of molecules of mrna for each technical replicate,  <dig> million random integers between  <dig> and z were sampled.
exon definitions
for this paper exons, rather than transcripts, were used.
estimating levels for individual transcript isoforms in rna-seq data can be quite complex as the observed data can be a mixture of isoforms  <cit> .
transcript assembly and abundance estimation from rna-seq reveals thousands of new transcripts and switching among isoforms.
to focus on the issue of technical replication, isoform estimation was not undertaken.
exons as defined by flybase  version  <dig>   <cit>  were mapped onto the genome.
overlapping exons in the same genome region were combined.
a single result per genomic region was recorded for each exon in that region .
for example, for exons with alternative start/end sites the longest genomic region was used to represent an exon.
the mapping between exons and the genome region was done via a perl script .
<dig>  of the  <dig>  genomic position exons represent uncomplicated annotation regions.
however, to ensure that results are not dependent on the exon definition, the observed data were used to derive exons.
consensus contigs were defined as contiguous regions with a per nucleotide read depth greater than or equal to  <dig> for all lanes within a particular experiment .
this definition results in the evaluation of genomic regions for which rna-seq data are available for all replicates of a particular experiment.
consensus contigs and their positions for these data are listed in additional files  <dig>   <dig>   <dig>
mapping reads
reads spanning junctions pose a problem.
each lane of data was first evaluated with tophat  <cit>  to create a set of read-supported junctions.
bowtie  <cit>  was used as the main tool for read alignment.
alignments were run with the "best" option to ensure that reported singleton alignments have the optimal score.
the "tryhard" flag was also used to induce bowtie to spend more time backtracking for hard-to-align reads.
reads aligning more than  <dig> times were removed from consideration.
first, a paired-end alignment was run.
unaligned reads from this step were aligned unpaired.
remaining unaligned reads were mapped to the junctions from tophat.
for d. simulans samples, reads were also aligned to d. simulans genomic reference sequence  <cit> .
bowtie is fast but cannot handle in/dels or reads that align partially to the reference.
reads remaining unaligned at this point were aligned with last  <cit> .
results were combined into a single alignment file in sam format and reported versus single reference .
analysis of read coverage was done with samtools  <cit> .
for each lane the resulting pileup file was the basis for all further calculations and the remaining unmapped reads were discarded.
detection and quantitation
each exon was projected onto each pileup file.
if at least one read maps to a particular exon then that exon was said to be detected.
rpkm which is the *c/mappable reads), where c is an estimate of the maximum number of reads  <cit>  was used to estimate the quantity of each exon.
analysis
to understand variability among replicates it is necessary to understand how well replicates agree, not how well they correlate  <cit> .
agreement between categorical variables can be calculated as the proportion of observations that share the same category divided by the total number of observations  <cit> .
if the number of observations in each category is unequal then two measurements may agree by chance more frequently than estimated using a simple agreement.
kappa calculates the agreement between two replicates and adjusts for chance agreement  <cit> .
the kappa statistic gives a measure of agreement and the value  <dig> represents perfect agreement  <cit> .
the most basic question is whether exons can be consistently detected between replicates.
the agreement of technical and biological replicates on detection of exons was calculated using both simple agreement and kappa statistics.
for ordinal scales, a weighted kappa allows for distance in categories to be used in the assessment of agreement and is interpreted in the same way as the kappa  <cit> .
the rpkm was grouped into  <dig> levels undetected/low-medium/high where high levels were defined as rpkm >  <dig>  alternatively rpkm was grouped into  <dig> levels using a log scale as follows: zero reads, rpkm less than  <dig>   <dig>   <dig>   <dig>   <dig>   <dig>   <dig> and greater than  <dig> were the ordinal categories used.
simple agreement  and weighted kappa statistics were calculated.
examination of the frequency table between the technical replicates can give insight into whether the disagreements are solely a function of the level of digital gene expression.
correlation coefficients have often been used to describe agreement.
however, correlation does not reflect agreement  <cit> .
spearman's correlation coefficient quantifies how well the relationship between two variables, , can be described using a monotonic function.
if one estimate is exactly twice the size of the other estimate then the correlation will be  <dig>  however, these two measures do not agree.
in addition to the scale factor, correlations depend upon the range of values observed.
a large range of values will produce a higher correlation than a narrow range of values.
in comparisons of whole genome expression the range is extremely large.
a further problem is the large number of data points makes the linear trend unlikely to be influenced by strong disagreement for some exons.
bland-altman  <cit>  proposed a method for visualizing data and providing insights into agreement.
these plots are known as mva or ma  plots in the gene expression literature and have been utilized intensively to assess presence of abundance-dependent bias  <cit> .
the plot is defined as follows: on the x axis is the average of the two measures  and on the y axis is the difference between the two measures.
these plots were constructed for exons for all combinations of technical replicates where there are non-zero measurements in both replicates.
an absolute difference of  <dig> has a different implication for agreement if the rpkm is  <dig> than if the rpkm is  <dig>  for that reason, bland-altman plots were also made on the natural log scale.
that is, the natural log of the rpkm was calculated for each technical replicate where the scale is considered to be additive rather than multiplicative and then used for computing the average and the difference.
coverage plots  <cit>  were constructed to examine whether exon length was associated with missing values.
to formally test whether exon abundance and/or length is related to whether an exon is present in all or only a portion of the technical replicates, an analysis of the dataset for the cell line c  <dig> was conducted.
only exons that were observed in at least one technical replicate were considered.
the dichotomous response of "in all technical replicates" versus "not in all technical replicates" was modeled using logistic regression.
formally, let yi =  <dig> if the ith exon is not present in all technical replicates,  <dig> otherwise.
the model
was fit, where ai is the abundance of the ith exon and li is the length of the ith exon and i =  <dig> ..,n where n is the number of exons observed in at least one technical replicate.
the average rpkm, taken across all technical replicates was used as the measure of abundance for that exon and the length of the exon in base pairs was used as the measure of exon length.
RESULTS
for the three experiments, the number of mappable reads per lane varies among the technical replicates and the number of exons detected increases with the number of mappable reads .
the exons that are detected vary between technical and biological replicates .
many exons are observed/present in one technical replicate but unobserved/absent in the other technical replicate.
there was no exception to this observation in any comparison made .
the agreement  in detection between technical replicates ranges from  <dig> - <dig> .
lanes with higher coverage have better agreement in detection, but even with the higher coverage rates, the number of exons detected varies.
the comparison with the least discrepancy had ~ <dig>  exons missing in one of the two technical replicates and several had more than  <dig>  missing .
among biological replicates discrepancies were larger than for technical replicates in that same experiment .
the simulation study shows that this variability in detection can be explained due to random noise from a very small sampling fraction .
a second simulation  based upon the multinomial model came to the same conclusion.
in the simulation, and in the real data examined here, discrepancies in detection are largely due to exons with coverage per nucleotide of less than  <dig> in the technical replicate in which they are detected .
this indicates that higher coverage increases the consistency of detection across technical replicates.
unfortunately, the vast majority of exons are observed at low coverage .
br indicates a biological replicate and tr indicates a technical replicate.
the experiments are described in figure  <dig>  there are a total number of  <dig>  exons corresponding to distinct genomic regions in flybase  <dig>  in this replicate  <dig>  replicates disagree by at least  <dig> log and  <dig> disagree by  <dig> or more logs.
all technical replicates were compared between biological replicates.
for example, br <dig> tr <dig> was compared to br <dig> tr <dig>  br <dig> tr <dig>  and br <dig> tr <dig> for all possible pairwise comparisons.
agreement in whether exons are detected , detected at a low level , and detected at a high level  and agreement on a  <dig> level ordinal scale as follows: rpkm less than  <dig>   <dig>   <dig>   <dig>   <dig>   <dig>   <dig> and greater than  <dig>  agreement was measured using a kappa coefficient.
the number of exons where disagreement is greater than  <dig> logs and greater than one log are also given.
agreement for common contigs is in additional file  <dig>  supplementary table s3
additionally, the length of the exon may contribute to the discrepancy.
to examine this, coverage plots were constructed .
coverage plots have a solid line if the exon was detected in the corresponding biological and technical replicate and a blank space if the exon was not present.
the plots were conducted by sorting the exons by length.
subsequent plots were made for exons surrounding a given percentile of exon length.
for exons that are detected inconsistently among replicates there seems to be no trend associated with the length.
plots of the 10th percentile appear to be similar to the 20th, 50th and 90th percentile .
the coverage plots do not take mrna abundance into account.
to test whether exon length and abundance, and the interaction between them contribute to differences in detection among technical replicates a logistic regression model was used.
the interaction between length and abundance was significant .
this means that the relationship between length and abundance differs between highly abundant and infrequently observed exons.
for infrequently observed exons, longer exons are more likely to be consistently detected than shorter exons.
this pattern is not observed in highly abundant exons.
it must be noted that this is a trend, and that some long exons are missing among some technical replicates and some short exons are present in all.
although main effects are often difficult to describe in the presence of interactions, in this case, it is clear that regardless of the exon length lower levels of abundance are less likely to be consistently detected .
within each lane, for each exon, the average number of reads per nucleotide  and the standard deviation were computed.
the coefficient of variation  versus the apn  was plotted for each exon.
coverage below apn  <dig> is highly variable.
the cv for an exon with a coverage of less than  <dig> reads per nucleotide is often greater than  <dig> and can be greater than  <dig>  the same pattern holds for all samples examined , indicating that this is not an artifact of the number of mappable reads.
in fact, when reads from the c <dig> experiment are combined and the same plot is constructed, the same pattern is seen .
although the maximum cv is smaller for the data combined across  <dig> lanes, than for a single lane of data, the trend is still apparent.
this, together with the results on consistency of detection, indicates that an average coverage of  <dig> reads is a useful target when planning an rna-seq study.
for a lane with almost  <dig> million mappable reads, not quite 60% of the exons in flybase were observed at this level of coverage.
there is great interest in understanding the performance of rna-seq and to that end many studies have used the correlation coefficient as a means of assessing concordance.
we plot the rpkms of technical replicates here against each other  to demonstrate that these plots are similar to those previously reported  <cit> .
as expected, pearson and spearman correlation coefficients are very high.
spearman correlation coefficients are in excess of  <dig> .
however, correlation is not a measure of agreement  <cit> .
by categorizing expression into several ordinal groups and then comparing two technical replicates , it becomes immediately obvious that there are areas of disagreement between the technical replicates, and these are not limited to exons with low coverage .
for one particular example given , the spearman correlation is  <dig>  and the weighted kappa statistic is  <dig> .
although the actual agreement will depend explicitly upon the ordinal values chosen and the scale used for expression, weighted kappa statistics for a range of possible categorizations were computed and found to be always lower than  <dig> , often close to  <dig>  and occasionally close to  <dig> .
this means that the use of rpkm as a measure of expression can be a log  different at all levels of expression.
this is not sensitive to the rpkm normalization, as normalization factors  <cit>  among technical replicates are close to  <dig>  indeed, this discrepancy is apparent in the simulated data, where normalization is not necessary.
sampling variance can result in differences in estimated expression of  <dig> log or greater.
this finding may explain some of the studies which report discordance between estimates of expression from rna-seq and known concentrations  <cit> .
rpkm was grouped into  <dig> categories on an approximate log <dig> scale.
the categories were: zero reads, average less than  <dig>   <dig>   <dig>   <dig>   <dig>   <dig>   <dig> and greater than  <dig>  technical replicates  <dig> and  <dig> were compared.
values on the diagonal are where the ordinal categories agree and off diagonal values  are disagreements.
the mva plot  <cit> , initially proposed by bland and altman  <cit>  to assess agreement between two methods of measuring an endpoint, is an intuitive examination of agreement that helps diagnose the magnitude and functional form of disagreements.
the bland-altman plot of the natural log-transformed data for the technical replicates of d. simulans biological replicate  <dig>  clearly shows that at lower levels of expression there is larger disagreement between technical replicates.
however, although the absolute disagreement is a function of abundance, the moving average smoother indicates that average disagreement or bias between technical replicates is consistently linear  over the abundance range.
this is in contrast to what is observed for microarray data where a nonlinear bias  is typically observed as a function of mean abundance making nonlinear normalization necessary  <cit> .
this is also consistent with the findings of bullard et al.
<cit>  who find that no nonlinear normalizations are needed.
does the choice of the flybase definition of an exon impact the results?
if the empirically defined common contig is used instead of the flybase exon, then by definition we are only examining regions seen in all technical replicates.
common contigs show greater levels of disagreement in the estimation of amount of expression between technical replicates than the flybase exons.
like the flybase exons, areas of lower coverage are more likely to disagree .
the disagreement across technical replicates decreases as the coverage increases for all the data examined here .
consistent with other reports  <cit>  biological variation is larger than technical variation.
the number of exons detected in both replicates is lower between biological replicates .
some of these differences may be due to the flow cell/lane assignments, which are unknown for all of these experiments.
consistent with previous reports, the level of coverage makes a large difference and some of the biggest disagreements are in the d. simulans experiments where  <dig> biological replicate has ~ <dig> million reads where the other biological replicates had ~ <dig> million reads.
these samples must have been run not only on different flow cells but also on different iterations of the technology.
based on previous findings, the finding that these samples have the largest differences and lowest agreement is not surprising.
discussion
the number of exons detected, was approximately 64%  in lanes with 5- <dig> million reads and 84% in lanes with approximately  <dig> million reads.
the results of this study indicate exons are not consistently detected among technical replicates when the average coverage for that exon is less than  <dig> reads per nucleotide.
additionally, the coefficient of variation can be very high for exons with coverage of less than  <dig> reads per nucleotide.
between lanes, exons that are detected in both lanes with one of those being at low coverage, disagree in the estimates of abundance  between technical replicates, leading us to conclude that when analyzing data for differential expression, one could consider identifying exons with coverage less than  <dig> for closer examination.
these findings are consistent across different versions of the technology, different species, and different laboratories .
the number of exons covered with an average of  <dig> reads per nucleotide--only 21% of all exons in lanes with 5- <dig> million reads and 58-60% in lanes with  <dig> million reads--is substantially smaller than the number of exons detected.
although agreement improves once coverage is above  <dig>  it is not perfect and disagreement can be large even when coverage is high.
why are there such large discrepancies in detection and in estimates of expression across technical replicates?
the underlying assumption for all gene expression technology, whether microarray, sage, q-pcr or rna-seq, is that the final measurement of gene expression for a particular sample is proportional to the underlying population of mrna extracted from that sample.
the variability in exon detection, and the disagreements in estimates of expression are perplexing.
a possible explanation would be the choice of the particular normalization.
however, the results are virtually identical if we use an average per nucleotide or rpkm.
a normalization constant  <cit>  was also estimated and found to be close to  <dig> among technical replicates, indicating that these results cannot be explained by improper normalization.
there are at least two obvious explanations for observing this disagreement: insufficient mixing and a very low sampling fraction.
if any of the solutions generated during preparation for cluster generation are insufficiently mixed, then the distribution of reads across samples is expected to be uneven.
this would be similar to what student saw in his experiment with yeast cells  <cit> .
in this experiment, insufficient mixing was partially attributed to the tendency of yeast cells to "stick together in groups which was not altogether abolished even by vigorous shaking."
student developed the negative binomial in response to this observation.
in the negative binomial distribution the variance is larger than in a poisson distribution.
the standard poisson has a variance equal to the mean.
poisson distributions that allow for over dispersion are general, and the negative binomial is one particular type of over dispersion.
the overdispersed poisson has been proposed for modeling rna-seq data  <cit> .
for both the d. melanogaster and d. simulans experiments, technical replicates were derived from an individual library analyzed on multiple lanes.
thus, discordance between the technical replicates must have occurred after library synthesis, presumably during generation of the templates on the flow cell.
we can envision other experimental plans in which technical replicates are derived from an individual sample, prior to library synthesis.
in this latter case, discordance between technical replicates could arise during or after library preparation.
modifications of the protocol have been developed that omit potential bias due to amplification  <cit> .
an inspection of the cluster generation protocol , carried out after library synthesis, reveals additional places where variation might be introduced.
the initial step of the protocol entails at least one large dilution step.
any uneven distribution of molecules in the library at this point, whether due to incomplete mixing or dna-dna interactions assisted by cations  <cit> , could lead to differences in the molecules loaded onto a lane.
single stranded libraries may have better mixing.
if there is aggregation of molecules in solution, then the current protocol could be modified to a dilution series with particular care taken to ensure molecules are mixed well and neither aggregate nor degrade in the process.
non-uniformity in solution is not limited to cdna given that rna molecules, used as input into library synthesis, can form intra- and intermolecular interactions.
rna solutions containing high concentrations of divalent metal ions are particularly prone to forming rna aggregates  <cit> .
as an extreme example, more than half of the rna in a pool of random rna sequences was found to aggregate together through self-complementarity and precipitate out of solution under conditions favoring intramolecular interactions  <cit> .
although it is not possible to determine from these experiments if unequal mixing is a problem, removing or reducing the potential for this problem should be a relatively straightforward modification of the current protocols.
the other possibility is that discrepancies among technical replicates are due to the sampling fraction,  <dig> % for  <dig> million reads .
this fraction is so low that the behavior of random sampling in this particular scenario should be examined.
the simulation study performed here indicates that data are consistent with this hypothesis.
marioni   <cit>  also reported data consistent with randomness.
agreement between technical replicates in the simulated data  shows the same patterns as the three sets of real data presented here.
in the simulated data, there are also examples of genes with substantial disagreement among technical replicates.
this indicates that the observed disagreement can be explained by the low sampling fraction, and that the results are consistent with the expectations of random variation.
what is initially surprising is the size of the random fluctuation.
however, when the sampling fraction is considered, the results are less surprising.
this study demonstrates that the current nature of the technology is to expect variation between single lanes both in exon detection and estimation of gene expression using rpkm.
technical variation, while smaller than the biological variation, cannot be ignored and should be accounted for in the study design.
the fundamental principles of statistical experimental design are thus just as important as ever, even with modern technology.
r.a. fisher defined these fundamentals to be randomization, replication and blocking  <cit> .
the optimal experimental design strategy will depend on the objectives of the study at hand.
randomization and blocking are important to incorporate in order to avoid confounding of biological with systematic experimental effects.
likely experimental effects in mrna seq data include lane, flow cell or library preparation batch.
randomized block designs have been shown to be an efficient use of resources while ensuring legitimate comparisons to be made.
multiplex designs  <cit>  represent another good solution which would ensure more complete coverage of exons across biological replicates and eliminate lane as a potential confounding variable in small experiments.
consider the following simple example of three biological replicates for treatments and three biological replicates for controls for a drosophila experiment.
if each sample is run on a single lane and more coverage is needed, then another full set of  <dig> lanes is necessary.
in contrast, a single multiplexed run where all  <dig> samples are run on all  <dig> lanes has been proposed  <cit> .
exons seen in some technical replicates but not others are clearly a result of technical variation.
in this design the total coverage remains almost the same as in the simple design.
the added benefit is that small numbers of extra lanes can be run if more coverage is desired.
a replication strategy must be chosen as well with an appropriate balance between biological replicates and technical replicates  <cit> .
while, increasing the number of biological replicates increases the precision and generalizability of a study more than increasing the number of technical replicates; for studies where low abundant mrnas are the focus, increasing technical replication may also be important.
CONCLUSIONS
rna-seq experiments need to replicate the results both technologically and biologically as the technical variation in exon presence and absence as well as amount of coverage is not negligible.
consistent with random variation due to low sampling fractions, technical variation is most pronounced in the variability of detection of exons when coverage is low.
however, disagreements between estimates of expression can occur at all levels of coverage.
authors' contributions
lmm designed the study, did analyses and wrote the paper, kkl did analyses, designed and executed simulation study, contributed to writing of the paper, amm designed experiments and writing of the paper, va carried out all bioinformatics pipeline design and execution, alo designed experiments and contributed to the paper, ljy designed experiments and contributed the paper, and svn designed experiments and contributed to writing of the paper.
all authors read and approved the final manuscript.
supplementary material
additional file 1
overlapping exons combined into single genomic region.
d. melanogaster ovo gene  used as an example of combining overlapping exons into a single genome region for mapping purposes.
format pdf.
view with adobe.
click here for file
additional file 2
chromosome postions for overlapping exons combined into single genomic regions.
format tsv.
view with wordpad.
click here for file
additional file 3
perl script for mapping between exons and a genome region.
format pl.
view with wordpad.
click here for file
additonal file 4
perl script that converts contiguous sequences in a samtools consensus pileup to fasta, bed, or coordinate tables.
format pl.
view with wordpad.
click here for file
additional file 5
perl script to find common sequences from contig beds.
format pl.
view with wordpad.
click here for file
additional file 6
st_2a_d.melanogaster_common_contigs.bed.
bed file containing chromosome postion, start, end and common_contig_id for d. melanogaster experiment .
format bed.
view with wordpad.
click here for file
additional file 7
st_2b_dsimulans_common_contigs.bed.
bed file containing chromosome postion, start, end and common_contig_id for d. simulans .
format bed.
view with wordpad.
click here for file
additional file 8
st_2c_167_common_contigs.bed.
bed file containing chromosome postion, start, end and common_contig_id for c <dig> experiment .
format bed.
view with wordpad.
click here for file
additional file 9
agreement between technical replicates and biological replicates for rpkm measured on common contigs.
a comparison in the agreement of the estimation of the amount of expression.
agreement in whether common contigs are expressed, expressed at a low level , and expressed at a high level  and agreement on a  <dig> level ordinal scale as follows: rpkm less than  <dig>   <dig>   <dig>   <dig>   <dig>   <dig>   <dig> and greater than  <dig>  agreement was measured using a kappa coefficient.
the number of common contigs where disagreement is greater than  <dig> logs and greater than one log are also given.
format xls.
view with excel.
click here for file
additional file 10
length of the exon does not explain disagreement in technical replicates.
coverage plots of the c <dig> cell lines data.
the y axis is the average coverage across all technical replicates.
a bar is drawn if the exon is present  in that technical replicate.
the 10th percentile represents that bottom 10% of the exons in length while the 90th percentile represents the top 10% of exons by length.
format pdf.
view with adobe
click here for file
additional file 11
bland altman plots for each biological replicate.
click here for file
acknowledgements
this work was supported by national institutes of health ; national science foundation .
michelle arbeitman and joe dunham are kindly acknowledged for the sharing of the d. melanogaster female head and d. simulans male head data.

 rna-seq has the potential to transform the way we study gene structure and expression. the ability to identify novel exons and splice sites is just the beginning. although there have been claims that rna-seq is more sensitive and has a larger dynamic range than a microarray, these claims are now being challenged. for all the optimism surrounding rna-seq, there is growing evidence that estimating gene expression in an rna-seq environment is not straightforward. in addition, technical variance, sequencing bias and mapping bias have been reported. the power and promise of rna-seq technology demand our attention. new papers on experimental design will lead the way to more thoughtful experimentation. one important component in careful study planning is an understanding of the different sources of variability so they can be accounted for in the experimental design. variability can occur at many levels, both biological and technical. biological variability is unaffected by the technology, as genotypes, individuals and even cell types vary regardless of how they are measured. bullard et al. examined the effects of library construction, flow cell and lane on detection of differential expression. the data used there was commercial grade rna and the samples compared expected to be very divergent. they concluded that the biological variation was larger than technical variation, underscoring the importance of including biological replicates in the study design. the observation that biological variability is large does not imply that technical variability should be ignored in the experimental design of a study. high levels of noise in a technology must be considered during experimental design. marioni et al. demonstrate clearly that technical variability is not different from what is expected due to random noise. yet, the magnitude of the sampling error is not explicitly discussed. one particular feature of deep sequencing is that even though there are millions of reads, the proportion of the mrna that is actually sequenced is low. given a solexa/illumina rna-seq library of nm concentration with a mean insert size of bp and a volume of ul, the number of molecules in a typical rna library is estimated to be × on the solexa/illumina gaiix, approximately million of the total possible molecules are sampled in a given lane. this represents approximately % of the total number of available molecules. even with the anticipated increase in the next generation of technology, the sampling fraction is will be less than %. the impact of such a low sampling fraction on the statistical properties of rna-seq is the focus of this manuscript. we address the following questions: 1) does a substantial amount of technical variability exist? 2) is the impact of technical variability the same for all levels of coverage? the study presented here considers technical variation in three experiments using solexa/illumina technology. although the number of exons consistently detected is improved by increasing the number of mappable reads, the detection of exons with low levels of coverage is inconsistent among technical replicates. inconsistent detection of exons is most pronounced for exons with average coverage of less than reads per nucleotide. in addition, there can be substantial disagreement in the estimated level of expression among technical replicates, even when the coverage is high. although the technical variability is not unexpected, technical variability cannot be ignored in rna-seq experiments. methods experiments experiment 1: three independent samples of d. melanogaster female heads were collected with each sample representing a unique pool of biological material. each sample was prepared according to manufacturer's instructions and then the same library was run on two lanes of a solexa/illumina flow cell, resulting in two technical replicates for each biological replicate, runs were base-pair paired end. experiment 2: three independent samples of d. simulans male heads were collected with each sample representing a unique pool of biological material. each sample was prepared according to manufacturer's instructions and then the same library was run on two lanes of a solexa/illumina flow cell, resulting in two technical replicates for each biological replicate, runs were base-pair paired end. experiment 3: one sample of d. melanogaster cell line c was run on lanes yielding technical replicates for biological replicate. data from these experiments were bp reads. the experimental design for experiments - are outlined in figure all samples were run using solexa/illumina paired end procedures. the relationship of lanes is unknown for all three experiments. the first two experiments were conducted by joe dunham and michele arbeitman at usc and sequenced at oregon state core facility while the third experiment was conducted by susan celniker at lawrence livermore national laboratory. simulation the number of molecules present in the library is estimated to be × this estimate was derived using the illumina protocol. specifically, the starting material was assumed to be ng of mrna that then resulted in ng of library with an insert size of 250bp in a volume of ul for a final library concentration of nm. the number of pmoles is therefore. the number of molecules is found by multiplying avogadro's number which is molecules/moles and adjusting for units to give × molecules. on the current solexa/illumina technology, the gaiix, approximately million of the total possible molecules are sampled in a given lane. this represents approximately % of the total number of available molecules for analysis. to understand how random sampling may produce variation in technical replicates in an rna-seq experiment, a small simulation of two technical replicates was conducted. we assume that the approximately × molecules of mrna were transcribed from genes. gene expression for the genes was modeled using a gamma density function with shape parameter and scale parameter this ensures that many genes will have a small number of molecules assigned to them and some genes will have a large number of molecules assigned to them. the range of the gamma density function was examined over the interval from to the integers from to represent bins of genes with similar relative levels of expression. to simulate relative expression for each of the i =..., individual genes, a random integer from a discrete uniform distribution between and were drawn and the value of the gamma density ti was calculated at si. the number of total molecules of mrna assigned to gene i, gi, was, ni= ti/Σ * × rounded to the nearest integer. Σ = z. due to the rounding, z was approximately, but not exactly equal to × for each gi with ni molecules a unique integer between and z was assigned, resulting in a unique association between an integer and a gene. to simulate the sampling of molecules of mrna for each technical replicate, million random integers between and z were sampled. exon definitions for this paper exons, rather than transcripts, were used. estimating levels for individual transcript isoforms in rna-seq data can be quite complex as the observed data can be a mixture of isoforms. transcript assembly and abundance estimation from rna-seq reveals thousands of new transcripts and switching among isoforms. to focus on the issue of technical replication, isoform estimation was not undertaken. exons as defined by flybase version were mapped onto the genome. overlapping exons in the same genome region were combined. a single result per genomic region was recorded for each exon in that region. for example, for exons with alternative start/end sites the longest genomic region was used to represent an exon. the mapping between exons and the genome region was done via a perl script. of the genomic position exons represent uncomplicated annotation regions. however, to ensure that results are not dependent on the exon definition, the observed data were used to derive exons. consensus contigs were defined as contiguous regions with a per nucleotide read depth greater than or equal to for all lanes within a particular experiment. this definition results in the evaluation of genomic regions for which rna-seq data are available for all replicates of a particular experiment. consensus contigs and their positions for these data are listed in additional files mapping reads reads spanning junctions pose a problem. each lane of data was first evaluated with tophat to create a set of read-supported junctions. bowtie was used as the main tool for read alignment. alignments were run with the "best" option to ensure that reported singleton alignments have the optimal score. the "tryhard" flag was also used to induce bowtie to spend more time backtracking for hard-to-align reads. reads aligning more than times were removed from consideration. first, a paired-end alignment was run. unaligned reads from this step were aligned unpaired. remaining unaligned reads were mapped to the junctions from tophat. for d. simulans samples, reads were also aligned to d. simulans genomic reference sequence. bowtie is fast but cannot handle in/dels or reads that align partially to the reference. reads remaining unaligned at this point were aligned with last. results were combined into a single alignment file in sam format and reported versus single reference. analysis of read coverage was done with samtools. for each lane the resulting pileup file was the basis for all further calculations and the remaining unmapped reads were discarded. detection and quantitation each exon was projected onto each pileup file. if at least one read maps to a particular exon then that exon was said to be detected. rpkm which is the *c/mappable reads), where c is an estimate of the maximum number of reads was used to estimate the quantity of each exon. analysis to understand variability among replicates it is necessary to understand how well replicates agree, not how well they correlate. agreement between categorical variables can be calculated as the proportion of observations that share the same category divided by the total number of observations. if the number of observations in each category is unequal then two measurements may agree by chance more frequently than estimated using a simple agreement. kappa calculates the agreement between two replicates and adjusts for chance agreement. the kappa statistic gives a measure of agreement and the value represents perfect agreement. the most basic question is whether exons can be consistently detected between replicates. the agreement of technical and biological replicates on detection of exons was calculated using both simple agreement and kappa statistics. for ordinal scales, a weighted kappa allows for distance in categories to be used in the assessment of agreement and is interpreted in the same way as the kappa. the rpkm was grouped into levels undetected/low-medium/high where high levels were defined as rpkm > alternatively rpkm was grouped into levels using a log scale as follows: zero reads, rpkm less than and greater than were the ordinal categories used. simple agreement and weighted kappa statistics were calculated. examination of the frequency table between the technical replicates can give insight into whether the disagreements are solely a function of the level of digital gene expression. correlation coefficients have often been used to describe agreement. however, correlation does not reflect agreement. spearman's correlation coefficient quantifies how well the relationship between two variables,, can be described using a monotonic function. if one estimate is exactly twice the size of the other estimate then the correlation will be however, these two measures do not agree. in addition to the scale factor, correlations depend upon the range of values observed. a large range of values will produce a higher correlation than a narrow range of values. in comparisons of whole genome expression the range is extremely large. a further problem is the large number of data points makes the linear trend unlikely to be influenced by strong disagreement for some exons. bland-altman proposed a method for visualizing data and providing insights into agreement. these plots are known as mva or ma plots in the gene expression literature and have been utilized intensively to assess presence of abundance-dependent bias. the plot is defined as follows: on the x axis is the average of the two measures and on the y axis is the difference between the two measures. these plots were constructed for exons for all combinations of technical replicates where there are non-zero measurements in both replicates. an absolute difference of has a different implication for agreement if the rpkm is than if the rpkm is for that reason, bland-altman plots were also made on the natural log scale. that is, the natural log of the rpkm was calculated for each technical replicate where the scale is considered to be additive rather than multiplicative and then used for computing the average and the difference. coverage plots were constructed to examine whether exon length was associated with missing values. to formally test whether exon abundance and/or length is related to whether an exon is present in all or only a portion of the technical replicates, an analysis of the dataset for the cell line c was conducted. only exons that were observed in at least one technical replicate were considered. the dichotomous response of "in all technical replicates" versus "not in all technical replicates" was modeled using logistic regression. formally, let yi = if the ith exon is not present in all technical replicates, otherwise. the model was fit, where ai is the abundance of the ith exon and li is the length of the ith exon and i =..,n where n is the number of exons observed in at least one technical replicate. the average rpkm, taken across all technical replicates was used as the measure of abundance for that exon and the length of the exon in base pairs was used as the measure of exon length.  for the three experiments, the number of mappable reads per lane varies among the technical replicates and the number of exons detected increases with the number of mappable reads. the exons that are detected vary between technical and biological replicates. many exons are observed/present in one technical replicate but unobserved/absent in the other technical replicate. there was no exception to this observation in any comparison made. the agreement in detection between technical replicates ranges from -. lanes with higher coverage have better agreement in detection, but even with the higher coverage rates, the number of exons detected varies. the comparison with the least discrepancy had ~ exons missing in one of the two technical replicates and several had more than missing. among biological replicates discrepancies were larger than for technical replicates in that same experiment. the simulation study shows that this variability in detection can be explained due to random noise from a very small sampling fraction. a second simulation based upon the multinomial model came to the same conclusion. in the simulation, and in the real data examined here, discrepancies in detection are largely due to exons with coverage per nucleotide of less than in the technical replicate in which they are detected. this indicates that higher coverage increases the consistency of detection across technical replicates. unfortunately, the vast majority of exons are observed at low coverage. br indicates a biological replicate and tr indicates a technical replicate. the experiments are described in figure there are a total number of exons corresponding to distinct genomic regions in flybase in this replicate replicates disagree by at least log and disagree by or more logs. all technical replicates were compared between biological replicates. for example, br tr was compared to br tr br tr and br tr for all possible pairwise comparisons. agreement in whether exons are detected, detected at a low level, and detected at a high level and agreement on a level ordinal scale as follows: rpkm less than and greater than agreement was measured using a kappa coefficient. the number of exons where disagreement is greater than logs and greater than one log are also given. agreement for common contigs is in additional file supplementary table s3 additionally, the length of the exon may contribute to the discrepancy. to examine this, coverage plots were constructed. coverage plots have a solid line if the exon was detected in the corresponding biological and technical replicate and a blank space if the exon was not present. the plots were conducted by sorting the exons by length. subsequent plots were made for exons surrounding a given percentile of exon length. for exons that are detected inconsistently among replicates there seems to be no trend associated with the length. plots of the 10th percentile appear to be similar to the 20th, 50th and 90th percentile. the coverage plots do not take mrna abundance into account. to test whether exon length and abundance, and the interaction between them contribute to differences in detection among technical replicates a logistic regression model was used. the interaction between length and abundance was significant. this means that the relationship between length and abundance differs between highly abundant and infrequently observed exons. for infrequently observed exons, longer exons are more likely to be consistently detected than shorter exons. this pattern is not observed in highly abundant exons. it must be noted that this is a trend, and that some long exons are missing among some technical replicates and some short exons are present in all. although main effects are often difficult to describe in the presence of interactions, in this case, it is clear that regardless of the exon length lower levels of abundance are less likely to be consistently detected. within each lane, for each exon, the average number of reads per nucleotide and the standard deviation were computed. the coefficient of variation versus the apn was plotted for each exon. coverage below apn is highly variable. the cv for an exon with a coverage of less than reads per nucleotide is often greater than and can be greater than the same pattern holds for all samples examined, indicating that this is not an artifact of the number of mappable reads. in fact, when reads from the c experiment are combined and the same plot is constructed, the same pattern is seen. although the maximum cv is smaller for the data combined across lanes, than for a single lane of data, the trend is still apparent. this, together with the results on consistency of detection, indicates that an average coverage of reads is a useful target when planning an rna-seq study. for a lane with almost million mappable reads, not quite 60% of the exons in flybase were observed at this level of coverage. there is great interest in understanding the performance of rna-seq and to that end many studies have used the correlation coefficient as a means of assessing concordance. we plot the rpkms of technical replicates here against each other to demonstrate that these plots are similar to those previously reported. as expected, pearson and spearman correlation coefficients are very high. spearman correlation coefficients are in excess of. however, correlation is not a measure of agreement. by categorizing expression into several ordinal groups and then comparing two technical replicates, it becomes immediately obvious that there are areas of disagreement between the technical replicates, and these are not limited to exons with low coverage. for one particular example given, the spearman correlation is and the weighted kappa statistic is. although the actual agreement will depend explicitly upon the ordinal values chosen and the scale used for expression, weighted kappa statistics for a range of possible categorizations were computed and found to be always lower than, often close to and occasionally close to. this means that the use of rpkm as a measure of expression can be a log different at all levels of expression. this is not sensitive to the rpkm normalization, as normalization factors among technical replicates are close to indeed, this discrepancy is apparent in the simulated data, where normalization is not necessary. sampling variance can result in differences in estimated expression of log or greater. this finding may explain some of the studies which report discordance between estimates of expression from rna-seq and known concentrations. rpkm was grouped into categories on an approximate log scale. the categories were: zero reads, average less than and greater than technical replicates and were compared. values on the diagonal are where the ordinal categories agree and off diagonal values are disagreements. the mva plot, initially proposed by bland and altman to assess agreement between two methods of measuring an endpoint, is an intuitive examination of agreement that helps diagnose the magnitude and functional form of disagreements. the bland-altman plot of the natural log-transformed data for the technical replicates of d. simulans biological replicate clearly shows that at lower levels of expression there is larger disagreement between technical replicates. however, although the absolute disagreement is a function of abundance, the moving average smoother indicates that average disagreement or bias between technical replicates is consistently linear over the abundance range. this is in contrast to what is observed for microarray data where a nonlinear bias is typically observed as a function of mean abundance making nonlinear normalization necessary. this is also consistent with the findings of bullard et al. who find that no nonlinear normalizations are needed. does the choice of the flybase definition of an exon impact the results? if the empirically defined common contig is used instead of the flybase exon, then by definition we are only examining regions seen in all technical replicates. common contigs show greater levels of disagreement in the estimation of amount of expression between technical replicates than the flybase exons. like the flybase exons, areas of lower coverage are more likely to disagree. the disagreement across technical replicates decreases as the coverage increases for all the data examined here. consistent with other reports biological variation is larger than technical variation. the number of exons detected in both replicates is lower between biological replicates. some of these differences may be due to the flow cell/lane assignments, which are unknown for all of these experiments. consistent with previous reports, the level of coverage makes a large difference and some of the biggest disagreements are in the d. simulans experiments where biological replicate has ~ million reads where the other biological replicates had ~ million reads. these samples must have been run not only on different flow cells but also on different iterations of the technology. based on previous findings, the finding that these samples have the largest differences and lowest agreement is not surprising. discussion the number of exons detected, was approximately 64% in lanes with 5- million reads and 84% in lanes with approximately million reads. the results of this study indicate exons are not consistently detected among technical replicates when the average coverage for that exon is less than reads per nucleotide. additionally, the coefficient of variation can be very high for exons with coverage of less than reads per nucleotide. between lanes, exons that are detected in both lanes with one of those being at low coverage, disagree in the estimates of abundance between technical replicates, leading us to conclude that when analyzing data for differential expression, one could consider identifying exons with coverage less than for closer examination. these findings are consistent across different versions of the technology, different species, and different laboratories. the number of exons covered with an average of reads per nucleotide--only 21% of all exons in lanes with 5- million reads and 58-60% in lanes with million reads--is substantially smaller than the number of exons detected. although agreement improves once coverage is above it is not perfect and disagreement can be large even when coverage is high. why are there such large discrepancies in detection and in estimates of expression across technical replicates? the underlying assumption for all gene expression technology, whether microarray, sage, q-pcr or rna-seq, is that the final measurement of gene expression for a particular sample is proportional to the underlying population of mrna extracted from that sample. the variability in exon detection, and the disagreements in estimates of expression are perplexing. a possible explanation would be the choice of the particular normalization. however, the results are virtually identical if we use an average per nucleotide or rpkm. a normalization constant was also estimated and found to be close to among technical replicates, indicating that these results cannot be explained by improper normalization. there are at least two obvious explanations for observing this disagreement: insufficient mixing and a very low sampling fraction. if any of the solutions generated during preparation for cluster generation are insufficiently mixed, then the distribution of reads across samples is expected to be uneven. this would be similar to what student saw in his experiment with yeast cells. in this experiment, insufficient mixing was partially attributed to the tendency of yeast cells to "stick together in groups which was not altogether abolished even by vigorous shaking." student developed the negative binomial in response to this observation. in the negative binomial distribution the variance is larger than in a poisson distribution. the standard poisson has a variance equal to the mean. poisson distributions that allow for over dispersion are general, and the negative binomial is one particular type of over dispersion. the overdispersed poisson has been proposed for modeling rna-seq data. for both the d. melanogaster and d. simulans experiments, technical replicates were derived from an individual library analyzed on multiple lanes. thus, discordance between the technical replicates must have occurred after library synthesis, presumably during generation of the templates on the flow cell. we can envision other experimental plans in which technical replicates are derived from an individual sample, prior to library synthesis. in this latter case, discordance between technical replicates could arise during or after library preparation. modifications of the protocol have been developed that omit potential bias due to amplification. an inspection of the cluster generation protocol, carried out after library synthesis, reveals additional places where variation might be introduced. the initial step of the protocol entails at least one large dilution step. any uneven distribution of molecules in the library at this point, whether due to incomplete mixing or dna-dna interactions assisted by cations, could lead to differences in the molecules loaded onto a lane. single stranded libraries may have better mixing. if there is aggregation of molecules in solution, then the current protocol could be modified to a dilution series with particular care taken to ensure molecules are mixed well and neither aggregate nor degrade in the process. non-uniformity in solution is not limited to cdna given that rna molecules, used as input into library synthesis, can form intra- and intermolecular interactions. rna solutions containing high concentrations of divalent metal ions are particularly prone to forming rna aggregates. as an extreme example, more than half of the rna in a pool of random rna sequences was found to aggregate together through self-complementarity and precipitate out of solution under conditions favoring intramolecular interactions. although it is not possible to determine from these experiments if unequal mixing is a problem, removing or reducing the potential for this problem should be a relatively straightforward modification of the current protocols. the other possibility is that discrepancies among technical replicates are due to the sampling fraction, % for million reads. this fraction is so low that the behavior of random sampling in this particular scenario should be examined. the simulation study performed here indicates that data are consistent with this hypothesis. marioni also reported data consistent with randomness. agreement between technical replicates in the simulated data shows the same patterns as the three sets of real data presented here. in the simulated data, there are also examples of genes with substantial disagreement among technical replicates. this indicates that the observed disagreement can be explained by the low sampling fraction, and that the results are consistent with the expectations of random variation. what is initially surprising is the size of the random fluctuation. however, when the sampling fraction is considered, the results are less surprising. this study demonstrates that the current nature of the technology is to expect variation between single lanes both in exon detection and estimation of gene expression using rpkm. technical variation, while smaller than the biological variation, cannot be ignored and should be accounted for in the study design. the fundamental principles of statistical experimental design are thus just as important as ever, even with modern technology. r.a. fisher defined these fundamentals to be randomization, replication and blocking. the optimal experimental design strategy will depend on the objectives of the study at hand. randomization and blocking are important to incorporate in order to avoid confounding of biological with systematic experimental effects. likely experimental effects in mrna seq data include lane, flow cell or library preparation batch. randomized block designs have been shown to be an efficient use of resources while ensuring legitimate comparisons to be made. multiplex designs represent another good solution which would ensure more complete coverage of exons across biological replicates and eliminate lane as a potential confounding variable in small experiments. consider the following simple example of three biological replicates for treatments and three biological replicates for controls for a drosophila experiment. if each sample is run on a single lane and more coverage is needed, then another full set of lanes is necessary. in contrast, a single multiplexed run where all samples are run on all lanes has been proposed. exons seen in some technical replicates but not others are clearly a result of technical variation. in this design the total coverage remains almost the same as in the simple design. the added benefit is that small numbers of extra lanes can be run if more coverage is desired. a replication strategy must be chosen as well with an appropriate balance between biological replicates and technical replicates. while, increasing the number of biological replicates increases the precision and generalizability of a study more than increasing the number of technical replicates; for studies where low abundant mrnas are the focus, increasing technical replication may also be important.  rna-seq experiments need to replicate the results both technologically and biologically as the technical variation in exon presence and absence as well as amount of coverage is not negligible. consistent with random variation due to low sampling fractions, technical variation is most pronounced in the variability of detection of exons when coverage is low. however, disagreements between estimates of expression can occur at all levels of coverage. authors' contributions lmm designed the study, did analyses and wrote the paper, kkl did analyses, designed and executed simulation study, contributed to writing of the paper, amm designed experiments and writing of the paper, va carried out all bioinformatics pipeline design and execution, alo designed experiments and contributed to the paper, ljy designed experiments and contributed the paper, and svn designed experiments and contributed to writing of the paper. all authors read and approved the final manuscript. supplementary material additional file 1 overlapping exons combined into single genomic region. d. melanogaster ovo gene used as an example of combining overlapping exons into a single genome region for mapping purposes. format pdf. view with adobe. click here for file additional file 2 chromosome postions for overlapping exons combined into single genomic regions. format tsv. view with wordpad. click here for file additional file 3 perl script for mapping between exons and a genome region. format pl. view with wordpad. click here for file additonal file 4 perl script that converts contiguous sequences in a samtools consensus pileup to fasta, bed, or coordinate tables. format pl. view with wordpad. click here for file additional file 5 perl script to find common sequences from contig beds. format pl. view with wordpad. click here for file additional file 6 st_2a_d.melanogaster_common_contigs.bed. bed file containing chromosome postion, start, end and common_contig_id for d. melanogaster experiment. format bed. view with wordpad. click here for file additional file 7 st_2b_dsimulans_common_contigs.bed. bed file containing chromosome postion, start, end and common_contig_id for d. simulans. format bed. view with wordpad. click here for file additional file 8 st_2c_167_common_contigs.bed. bed file containing chromosome postion, start, end and common_contig_id for c experiment. format bed. view with wordpad. click here for file additional file 9 agreement between technical replicates and biological replicates for rpkm measured on common contigs. a comparison in the agreement of the estimation of the amount of expression. agreement in whether common contigs are expressed, expressed at a low level, and expressed at a high level and agreement on a level ordinal scale as follows: rpkm less than and greater than agreement was measured using a kappa coefficient. the number of common contigs where disagreement is greater than logs and greater than one log are also given. format xls. view with excel. click here for file additional file 10 length of the exon does not explain disagreement in technical replicates. coverage plots of the c cell lines data. the y axis is the average coverage across all technical replicates. a bar is drawn if the exon is present in that technical replicate. the 10th percentile represents that bottom 10% of the exons in length while the 90th percentile represents the top 10% of exons by length. format pdf. view with adobe click here for file additional file 11 bland altman plots for each biological replicate. click here for file acknowledgements this work was supported by national institutes of health; national science foundation. michelle arbeitman and joe dunham are kindly acknowledged for the sharing of the d. melanogaster female head and d. simulans male head data.
technical variability is too high to ignore. technical variability results in inconsistent detection of exons at low levels of coverage. we provide practical recommendations for dealing with the technical variability, without dramatic cost increases. biological variation has been reported to be larger than technical variation. the size of technical variance, and the role of sampling are examined in this manuscript. in addition, technical variation has been reported to be in line with expectations due to random sampling. however, strategies for dealing with technical variation will differ depending on the magnitude. exon detection between technical replicates is highly variable when the coverage is less than reads per nucleotide and estimates of gene expression are more likely to disagree when coverage is low. in this study three independent solexa/illumina experiments containing technical replicates are analyzed. when coverage is low, large disagreements between technical replicates are apparent. this may be due to the low sampling fraction and if so, it will persist as an issue needing to be addressed in experimental design even as the next wave of technology produces larger numbers of reads.
 rna-seq is revolutionizing the way we study transcriptomes. mrna can be surveyed without prior knowledge of gene transcripts. alternative splicing of transcript isoforms and the identification of previously unknown exons are being reported. initial reports of differences in exon usage, and splicing between samples as well as quantitative differences among samples are beginning to surface. biological variation has been reported to be larger than technical variation. in addition, technical variation has been reported to be in line with expectations due to random sampling. however, strategies for dealing with technical variation will differ depending on the magnitude. the size of technical variance, and the role of sampling are examined in this manuscript.  in this study three independent solexa/illumina experiments containing technical replicates are analyzed. when coverage is low, large disagreements between technical replicates are apparent. exon detection between technical replicates is highly variable when the coverage is less than reads per nucleotide and estimates of gene expression are more likely to disagree when coverage is low. although large disagreements in the estimates of expression are observed at all levels of coverage.  technical variability is too high to ignore. technical variability results in inconsistent detection of exons at low levels of coverage. further, the estimate of the relative abundance of a transcript can substantially disagree, even when coverage levels are high. this may be due to the low sampling fraction and if so, it will persist as an issue needing to be addressed in experimental design even as the next wave of technology produces larger numbers of reads. we provide practical recommendations for dealing with the technical variability, without dramatic cost increases. 