BACKGROUND
analysis of the expression of mrna produces a molecular readout of the biological activity in a tissue or entire organism sample, reflecting which parts of the genome are being transcribed and how much of each transcript is available for translation into protein or to perform a regulatory role.
comparing mrna expression across different conditions, such as developmental stages or after gene perturbation, helps to unravel the complexities of biological systems.
in recent years high-throughput short tag rna sequencing technology  has provided a convenient tool for delving deeper into mrna expression  <cit>  using the whole or part of each transcript .
more specifically, sequence reads are converted into count data with the aim of quantifying transcriptional differences between biological samples using polya pull down at the 3′ end of transcripts  and cap analysis gene expression  at the 5′ end of transcripts  <cit> .
many rna-seq methods are excellent tools for in-depth mrna expression analysis of small numbers of samples and provide information on the entire rna molecule, alternative splicing and the quantity of transcript.
however, they require an involved library preparation and often complex sequence analysis  <cit>  and are not amenable to large-scale application with a fast turn-around.
we present a purely quantitative digital gene expression sample processing and analysis package called differential expression transcript counting technique  that begins with tissue samples and produces a text table or html table, comprising genomic coordinates representing the 3′ ends of genes, raw and normalised counts, and a fold change in transcript abundance between two conditions with an associated p-value.
our simplified library preparation and analysis protocol incorporates a sample indexing system and allows processing and sequencing of large numbers of samples and replicates.
the genomic coordinates can be compared to existing gene annotation, but they also identify unannotated genomic regions showing an alteration in polya+ transcript number.
to assess the utility of the pipeline we used zebrafish mutants carrying loss of function alleles from the zebrafish mutation project   <cit>  and compared morphologically abnormal embryos with normal sibling embryos.
RESULTS
library preparation and sequence processing
we selected four mutant zebrafish lines from the zebrafish mutation project  to test the differential expression transcript counting technique  pipeline.
we collected morphologically normal and abnormal sibling single embryos in replicates from the same clutch obtained from timed single pair matings to synchronize the developmental stage.
total rna was extracted from single zebrafish embryos with sufficient residual dna to confirm the genotype of each embryo by kasp genotyping  <cit> .
libraries were prepared from  <dig> ng of total rna.
several features make our libraries different to standard rna-seq methods  <cit> .
we have simplified the library preparation by reducing the number of steps, but have added several useful modifications .
the dnasei digestion has been combined with the rna fragmentation step and is followed by the first anchored polya pulldown enrichment.
while the rna molecules are immobilised on magnetic beads rna to rna ligation introduces part of illumina adapter sequence  <dig>  after elution we perform a second round of 3′ end enrichment with an anchored oligo dt reverse transcription primer.
this primer also incorporates a sample-specific in-read index sequence, a unique molecular identifier  sequence and part of illumina adapter  <dig>  after reverse transcription through the captured rna molecules and the partial illumina adapter  <dig> sequence, illumina adapters are completed during a final library amplification step.
the replicate libraries for each allele were pooled and sequenced on one lane equivalent  or two lane equivalents  by illumina hiseq  <dig>  figure 2a shows a sequence depth of between  <dig> and  <dig> million read pairs per lane equivalent, with  <dig> % and  <dig> % of the sequence mapping to the zv <dig> zebrafish reference genome, respectively.
figure 2a also shows a relatively even quantity of sequence per library with the occasional outlier.
during library preparation the sample amplification can result in duplicate reads for the same original transcript which is particularly relevant if an unknown amount of rna is accidentally lost before amplification.
the duplicate rate in a library reflects library complexity and is therefore an important quality indicator for each library.
duplicate reads caused by amplification can be identified by incorporating random sequences as unique molecular identifiers  into the primary sample .
we use a modified version of picard markduplicates  <cit>  and flag reads as potential duplicates if they share outer coordinates with other mapped read pairs and have the same umi.
figure 2b shows that accounting for the umis reduced the median duplicate rate from  <dig>  % to  <dig>  % with a few libraries showing a higher duplicate rate of up to  <dig> %.
if the duplicate rate goes above  <dig> % then we examine laboratory procedures for technical issues such as rnase contamination in reagents.
this method remains an estimate of library duplication due to the possibility of two independent molecules carrying the same umi and umi sequence alteration during subsequent amplification.
read  <dig> is used to predict a single genomic position defining the transcript counting 3′ end .
between  <dig> % and  <dig> % of the read 2s map to the reference genome and where they accumulate peaks are called and reads are quantified as counts .
read pair information attaches these count data to the tc 3′ end and the in-read index sequence identifies their sample origin.
deseq <dig>  <cit>  is used to estimate differential transcript abundance between conditions, independent of gene annotation.
while the strand-specific tc 3′ ends can be linked to any gene annotation, here we use the ensembl gene build  <cit> .
fig.
2d shows the total number of regions called as peaks , the subset associated with a gene where deseq <dig> has estimated an adjusted p-value , genes showing differential transcript abundance between conditions with an adjusted p-value <=  <dig>  and finally the subset of the latter showing differential transcript abundance with a fold change >  <dig>  the analysis pipeline utilises a single configuration file describing the samples, the location of the sequence files and the conditions, and with one command converts duplicate-marked bam files into the detct pipeline output gene list tables.
these tables show the closest strand-specific ensembl gene to the tc 3′ end, the region called as a peak, the unadjusted p-value, the adjusted p-value, the fold change between conditions and the count data.
it is presented in tab or comma-separated tables or as an html table .fig.
<dig> detct pipeline workflow.
between nine and  <dig> pairs of mutant and normal zebrafish embryos were collected from one clutch and rna extracted.
a following dnasei treatment and chemical fragmentation, molecules representing the 3′ end of transcripts were enriched by pulldown using an anchored biotinylated oligo dt primer attached to streptavidin magnetic beads .
an rna oligo matching part of the illumina read  <dig> adapter  was ligated onto the 5′ end, the rna eluted and annealed to an oligo comprising partial read  <dig> illumina adapter  followed by  <dig> random bases , then an eight base indexing sequence specific to each sample  and finally a  <dig> base anchored polyt sequence .
after reverse transcription the illumina adapter sequences were completed during library amplification.
libraries were quantified, pooled in equimolar amounts and sequenced by illumina hiseq  <dig>  b after decoding the indexing sequence, the trimmed zebrafish sequences  were mapped to the reference genome and duplicate reads were flagged.
c the coordinate representing the transcript counting 3′ end  was predicted using the base immediately 3′ of the polyt sequence in read  <dig> .
after calling peaks using all mapped read 2s the resulting counts were associated with their respective sample .
the count data were used to identify differential transcript abundance between conditions using deseq <dig>  <cit>  and reported as a fold change with an adjusted p-value.
the tc 3′ ends were matched to the closest ensembl transcript 3′ ends on the same strand .
gene list tables were produced and ordered by the lowest adjusted p-value.
these gene lists were filtered for genes showing differential transcript abundance using the adjusted p-value and the proximity of the tc 3′ end and ensembl gene end
fig.
<dig> library features.
four sets of libraries were prepared from zebrafish embryo lines carrying a disrupted gene.
morphologically abnormal and normal sibling embryos were collected from in-crossed lines and total rna extracted.
for each mutant line  <dig> to  <dig> samples of indexed libraries were made and sequenced by illumina hiseq  <dig>  a the number of reads and the number mapping to the zv <dig> reference genome per library are shown.
the total reads, the mean per library and the standard deviation are shown on the right.
b for each library  the proportion of reads identified as duplicates using outer mapping coordinates alone are shown on the x-axis and after accounting for the unique molecule identifier  on the y-axis.
c the reads were passed through the detct analysis pipeline.
the number of read 2s mapped and the number of counts called under peaks as discrete regions are shown per library sample.
the total reads, the mean per library and the standard deviation are shown on the right.
d using the gene list output from the detct pipeline the chart shows the number of discrete regions identified per collection of libraries, the number of genes with an adjusted p-value  from deseq <dig> , the number of genes with an adjusted p-value <=  <dig>  and the number of genes with an adjusted p-value <=  <dig>  plus a fold change  >  <dig>  these analyses were performed with a stringent  or relaxed  proximity filter between the ensembl transcript and tc 3′ end
data interpretation
the experimental rationale suggests that all transcript counting 3′ ends  should match the 3′ ends of transcripts.
however, gene annotation is sometimes incomplete and occasionally both the annotation and tc 3′ end can arise from experimental artefact.
therefore the detct output specifies the distance between the tc 3′ end and the nearest ensembl transcript 3′ end on the same strand.
naturally, this assumes the transcript to be annotated correctly and this is not necessarily true .
the ensembl transcript ends that exactly match the tc 3′ ends can be easily filtered from the results table.
non-exact matches suggest incomplete gene annotation or novel alternative transcript ends and both situations can be validated by individual inspection .
choosing a close proximity filter of the coordinates, such as between - <dig> bases  and + <dig> bases , reduces the likelihood of a false positive match and using these criteria we were able to detect a mean of  <dig> genes per experiment .
in contrast a more relaxed proximity filter, such as between - <dig> to + <dig>  identifies many more genes , but also finds more false positive ends.
this is discussed further in the comparison to rna-seq below.
one cause of false positive tc 3′ ends is oligo dt priming from polya or degenerate polya sequence within rna molecules or possibly from residual dna.
we identified falsely primed tc 3′ ends during the detct analysis pipeline by examining the  <dig> bases 3′ of the tc 3′ end and removed those which potentially derived from non-polya tail priming using the criteria described in the methods .
one zebrafish-specific example is in the mitochondrial genome where a region rich in adenine  in the rrna ensdarg <dig> accounts for between  <dig> % to  <dig> % of the total counts in the four test experiments .
this example escaped our current polya filtering method.
to mitigate false positive 3′ end calls, we have begun to prepare a list of tc 3′ ends we believe not to be true transcript tc 3′ ends.
similarly, we are able to catalogue all the true positive tc 3′ ends we find and build cross-experiment profiles of regularly identified tc 3′ ends.
we further filtered the results file by restricting the regions used to those with the 3′ coordinate of the region not more than  <dig> bases upstream of the tc 3′ end.
assessing the detct method
to assess the variation between libraries we extracted rna from a pool of zebrafish embryos and made  <dig> replicate transcript counting libraries with different indexing sequences using  <dig> μg each.
libraries were pooled, sequenced on an illumina miseq and analysed using the detct pipeline.
the number of normalised counts for each genomic region called by detct for each library was determined.
the regions were filtered for a maximum  <dig> bp distance to an ensembl gene 3′ end, as described in additional file  <dig>  and compared using a pearson correlation .
these  <dig> libraries show our method displays good technical reproducibility .fig.
<dig> technical replicate.
twelve replicate transcript counting libraries were prepared from  <dig> μg of a pool of wild-type zebrafish embryo total rna sample.
libraries were sequenced by illumina miseq and analysed using the detct pipeline.
a the normalised counts for each region were extracted .
after filtering the data for only the regions which we would use to call differential transcript abundance the counts from all  <dig> libraries were compared using a pearson correlation .
cells coloured yellow in the pearson correlation are the most highly correlated while those in blue are the least correlated with a colour gradient inbetween.
b in addition four concentrations of ercc spike mix  <dig> were added in triplicate to the same  <dig> libraries prior to library construction.
we added the quantity suggested by the manufacturer , five times the quantity , one fifth of the quantity  and one tenth of the quantity .
the reads were mapped to the zebrafish reference sequence and ercc spike reference sequence.
the diagram shows the  <dig> ercc spikes represented by a circle in descending order of spike copy number in the mix on the x-axis and spike abundance on the y-axis.
the blue circles show spikes identified in the detct pipeline while those in red were not found.
c the detct analysis pipeline was run using six libraries at a time with three replicates as one pair of conditions and in all six possible condition combinations.
the mean log <dig> fold change was calculated for all the spikes detected by the detct analysis and plotted against the expected log <dig> fold change as circles.
each circle is labelled with the conditions being compared and observed log <dig> fold change over the expected log <dig> fold change.
the numbers in brackets indicate how many spikes show differential transcript abundance
to examine the method’s performance in quantifying differential transcript abundance we added ercc spike mix  <dig>  to the total rna prior to making the  <dig> technical replicate libraries described above.
we added the quantity of spike mix recommended by the manufacturer, five times the quantity, a fifth of the quantity and a tenth of the quantity, all in triplicate, to create four conditions.
after analysis through the detct pipeline we identified  <dig> spikes .
the spikes with high copy number in the mix were detected, but not those with lower copy number .
the mean log <dig> fold changes were calculated for the spikes detected in all combinations of condition and compared to the log <dig> fold change expected .
we found a good correlation between observed and expected log <dig> fold change.
additionally, in three of the pairwise comparisons no false positives were detected and in the other three only a total of six non-ercc regions were found to be differentially expressed, suggesting that the method has high specificity.
our method was designed to maximise the number of tissue samples we can process with relatively shallow sequencing whilst still obtaining sufficient information to implicate gene networks modified by the condition change.
the ability to make numerous transcript counting libraries has two main advantages.
firstly, problem libraries resulting from sample loss or showing low complexity can be excluded from analysis and their removal has little effect on the statistical power of the analysis.
similarly, samples with incorrectly assigned condition have less influence on the final result.
secondly, increasing the number of replicates improves the statistical power of the analysis.
to assess the impact of increasing the number of replicates whilst retaining the same total amount of sequencing we performed a permutation test using the zmp_ph <dig> data  <cit> .
reads were combined to produce two samples comprising  <dig> , <dig> mutant and  <dig> , <dig> wild-type reads.
these were randomly split ten times into ten collections of pseudo-samples each containing equal numbers of reads but with the number of pairs of pseudo-samples ranging from two pairs to  <dig> pairs .
after passing these collections of pseudo-samples through the detct pipeline  we identified  <dig>  transcript counting 3′ ends  in all  <dig> simulations.
although this simulation has removed the variance from the original biological replicates it shows that increasing the number of libraries at the expense of read depth improves detection of tc 3′ ends showing differential abundance, as was previously noted  <cit> .fig.
<dig> sample number.
the  <dig> mutant and  <dig> wild-type zmp_ph <dig> bam files were downsampled and merged to create two bam files, one from the mutant samples and one from the wild-type samples.
the real index sequences in each bam file were then replaced with fake index sequences in order to simulate assigning the reads at random to between  <dig> to  <dig> different pseudo-samples.
each of these  <dig> permutations was repeated  <dig> times, to give  <dig> simulations in total , each of which was run through the detct pipeline.
the chart shows the average total number of counts for each pseudo-sample  and the number of  tc 3′ ends with a relaxed ensembl transcript proximity filter showing differential transcript abundance  with the mean above each group of ten
many mrna expression pipelines use whole transcript rna-seq protocols and a range of analysis tools  <cit> .
we don’t intend to replace these methods but present an alternative method for high-throughput mrna expression screening.
however, it is useful to compare the results of both methods.
to this end we made detct libraries and standard non-directional polya pulldown illumina rna-seq libraries for two alleles.
the same three wild-type and three mutant zebrafish total rna samples were processed for each method, plus two or six additional libraries using the detct protocol.
we sequenced one hiseq  <dig> lane equivalent for each allele by each method .
potential duplicate reads were identified and eliminated.
the umi in detct allowed a more accurate identification of duplicates and hence fewer reads were dismissed as duplicates.
read 2s from the rna-seq data that mapped to the genome were compared to ensembl gene annotation to produce count data for each gene and the detct pipeline was used to extract count data linked to ensembl transcripts from the transcript counting reads.
deseq <dig> was run on both sets of count data.
even after removing duplicate reads there are generally more counts in the rna-seq count data .
this is probably due to a drop in read quality following the oligo dt sequence in transcript counting  read  <dig>  rna-seq initially identifies more genes, however, the gap between the methods is reduced substantially when regions with a low mean count  are filtered out by deseq <dig> .
the number of tc 3′ ends  is higher than that of detct genes  which suggests alternative 3′ ends  <cit>  .
however, some represent false positive tc 3′ ends which escape our filter for false oligo dt priming.
to assess the genes showing differential transcript abundance the full gene list was filtered for protein-coding genes with an adjusted p-value of <=  <dig>  and an absolute fold change >=  <dig>  between mutant and wild-type .
rna-seq identifies more genes, but using a less stringent proximity filter between the tc 3′ end and the 3′ end of an ensembl transcript increases the detection rate in detct.
for a direct comparison between the two methods we identified the genes with an adjusted p-value <=  <dig>  and absolute fold change >=  <dig>  for both methods and applied the less stringent proximity filter to the tc 3′ ends.
the fold change for these genes in rna-seq and tc was compared  and showed a good correlation with r2 =  <dig>  , which suggests the two methods are finding the same alterations in transcript abundance.
we next looked at genes which show an adjusted p-value <=  <dig>  and an absolute fold change >=  <dig> by one method, but failed to meet one or both criteria by the other .
for both methods  <dig> genes have an absolute fold change >=  <dig>  but fail to have sufficient power to call an adjusted p-value <=  <dig>  by one method.
similarly, for  <dig> genes where one method fails to show an absolute fold change >=  <dig> the actual fold change is just below  <dig>  suggesting further genes where the two methods give comparable fold change results.
we then applied the stringent tc 3′ end proximity filter to the same data which led to the removal of  <dig> genes .
examining the tc 3′ end of these  <dig> genes showed they fell into two groups, either true tc 3′ ends or false tc 3′ ends assumed to be derived from experimental artefact .
where both methods gave an adjusted p-value <=  <dig>  and an absolute fold change >=  <dig> all  <dig> were shown to be true ends .
by contrast, in the gene sets only called by one method 11/ <dig> tc 3′ ends lost to more stringent filtering were false positives.
together this analysis shows the removal of  <dig> genes by increasing the stringency of the detct proximity filter resulted in losing  <dig> true positives , but prevented calling  <dig> false positive tc 3′ ends.table  <dig> comparison of rna-seq and transcript counting
lamc1
sa379
mdn1
sa1349
fig.
<dig> comparing rna-seq and transcript counting.
the genes identified as showing differential transcript abundance in experiment lamc1
sa <dig> with an adjusted p-value <=  <dig>  and an absolute fold change of >=  <dig> by one or both methods at the less stringent proximity filter between the transcript counting 3′ end  and the ensembl transcript end  were plotted with their log <dig> fold change from rna-seq against the log <dig> fold change from tc.
genes with an adjusted p-value <=  <dig>  and absolute fold change >=  <dig> in both rna-seq and tc are shown as blue circles.
genes which fail on one or both criteria in rna-seq but not in tc are shown as green circles and vice versa as red circles.
a the log <dig> fold changes of genes with the tc 3′ end and ensembl transcript proximity filter between - <dig> and + <dig>  the arrows highlight examples of genes which are not seen in graph b. b the log <dig> fold changes of genes with the proximity filter between - <dig> and + <dig>  c a table showing the number of genes represented by each circle colour.
the fourth and fifth columns show the genes which are lost during the more stringent proximity filtering and whether these are considered true positives or false positives after examining the tc 3′ end.
note that for two genes where the rna-seq and tc data correlated, the closest tc 3′ end was found to be false, but a true end was identified further downstream and indicated by square brackets.
the number in curly brackets indicates genes with a true tc 3′ end and an absolute fold change > <dig> but which fail to show an adjusted p-value <=  <dig>  in tc
CONCLUSIONS
we present a quantitative mrna transcript profiling package that starts with tissue samples and produces a gene list by counting the 3′ end of any polyadenylated transcripts using illumina sequencing.
unlike whole transcript rna-seq each transcript is counted only once giving a more representative estimate of transcript abundance  <cit> .
short, rare transcripts are as likely to be represented in the sequence as long, rare transcripts.
assaying only the 3′ end of transcripts is also more resilient to degraded rna samples, particularly if cells or tissues are compromised by a treatment.
differential transcript abundance is identified based on genome sequence and is independent of gene annotation.
this highlights regions containing novel transcripts as well as previously undescribed alternative transcript 3′ ends, which are implicated in biological processes related to the condition.
each alternative transcript 3′ end is represented by discrete count data and has the potential to add layers of functional annotation to sequence at the 3′ ends of gene models.
our streamlined library preparation reduces material loss and allows us to produce large numbers of libraries.
working with more replicates not only provides more power to the differential expression calculation but also lessens the impact of occasional sample loss or failure.
this allows us to measure transcriptional changes even in rare tissues within the whole organism.
the addition of a unique molecular identifier  helps us to assess the quality of each library preparation and we can remove any which underperform from the analysis.
simultaneously processing a large number of samples allows us to ask complex questions, such as what is the impact of different concentrations of a compound at varying stages of embryonic development?
similarly, we could screen large numbers of mutagenised individuals following an infection challenge to identify variation in infection response and inflammation.
in both cases relatively shallow sequencing of each sample will result in reduced sensitivity, but can highlight a critical combination of conditions to explore further using deeper sequencing of selected libraries to increase the number of transcripts detected or by an alternative method such as rna-seq.
we are also able to assay large numbers of wild-type samples producing a gene profile baseline that can be refined with data from each individual biological replicate.
for precious or difficult to obtain rna samples conventional full transcript rna-seq methods are more suited, but for rapid assessment of many samples with a desire to implicate a molecular process, transcript counting is the better option.
the gene lists provided by detct are the link between a living organism and the array of gene ontology, gene expression and gene interaction network data in the public domain.
we believe the power of the differential expression transcript counting technique  lies in effectively and efficiently bridging this gap.

 analysis of the expression of mrna produces a molecular readout of the biological activity in a tissue or entire organism sample, reflecting which parts of the genome are being transcribed and how much of each transcript is available for translation into protein or to perform a regulatory role. comparing mrna expression across different conditions, such as developmental stages or after gene perturbation, helps to unravel the complexities of biological systems. in recent years high-throughput short tag rna sequencing technology has provided a convenient tool for delving deeper into mrna expression using the whole or part of each transcript. more specifically, sequence reads are converted into count data with the aim of quantifying transcriptional differences between biological samples using polya pull down at the 3′ end of transcripts and cap analysis gene expression at the 5′ end of transcripts. many rna-seq methods are excellent tools for in-depth mrna expression analysis of small numbers of samples and provide information on the entire rna molecule, alternative splicing and the quantity of transcript. however, they require an involved library preparation and often complex sequence analysis and are not amenable to large-scale application with a fast turn-around. we present a purely quantitative digital gene expression sample processing and analysis package called differential expression transcript counting technique that begins with tissue samples and produces a text table or html table, comprising genomic coordinates representing the 3′ ends of genes, raw and normalised counts, and a fold change in transcript abundance between two conditions with an associated p-value. our simplified library preparation and analysis protocol incorporates a sample indexing system and allows processing and sequencing of large numbers of samples and replicates. the genomic coordinates can be compared to existing gene annotation, but they also identify unannotated genomic regions showing an alteration in polya+ transcript number. to assess the utility of the pipeline we used zebrafish mutants carrying loss of function alleles from the zebrafish mutation project and compared morphologically abnormal embryos with normal sibling embryos.  library preparation and sequence processing we selected four mutant zebrafish lines from the zebrafish mutation project to test the differential expression transcript counting technique pipeline. we collected morphologically normal and abnormal sibling single embryos in replicates from the same clutch obtained from timed single pair matings to synchronize the developmental stage. total rna was extracted from single zebrafish embryos with sufficient residual dna to confirm the genotype of each embryo by kasp genotyping. libraries were prepared from ng of total rna. several features make our libraries different to standard rna-seq methods. we have simplified the library preparation by reducing the number of steps, but have added several useful modifications. the dnasei digestion has been combined with the rna fragmentation step and is followed by the first anchored polya pulldown enrichment. while the rna molecules are immobilised on magnetic beads rna to rna ligation introduces part of illumina adapter sequence after elution we perform a second round of 3′ end enrichment with an anchored oligo dt reverse transcription primer. this primer also incorporates a sample-specific in-read index sequence, a unique molecular identifier sequence and part of illumina adapter after reverse transcription through the captured rna molecules and the partial illumina adapter sequence, illumina adapters are completed during a final library amplification step. the replicate libraries for each allele were pooled and sequenced on one lane equivalent or two lane equivalents by illumina hiseq figure 2a shows a sequence depth of between and million read pairs per lane equivalent, with % and % of the sequence mapping to the zv zebrafish reference genome, respectively. figure 2a also shows a relatively even quantity of sequence per library with the occasional outlier. during library preparation the sample amplification can result in duplicate reads for the same original transcript which is particularly relevant if an unknown amount of rna is accidentally lost before amplification. the duplicate rate in a library reflects library complexity and is therefore an important quality indicator for each library. duplicate reads caused by amplification can be identified by incorporating random sequences as unique molecular identifiers into the primary sample. we use a modified version of picard markduplicates and flag reads as potential duplicates if they share outer coordinates with other mapped read pairs and have the same umi. figure 2b shows that accounting for the umis reduced the median duplicate rate from % to % with a few libraries showing a higher duplicate rate of up to %. if the duplicate rate goes above % then we examine laboratory procedures for technical issues such as rnase contamination in reagents. this method remains an estimate of library duplication due to the possibility of two independent molecules carrying the same umi and umi sequence alteration during subsequent amplification. read is used to predict a single genomic position defining the transcript counting 3′ end. between % and % of the read 2s map to the reference genome and where they accumulate peaks are called and reads are quantified as counts. read pair information attaches these count data to the tc 3′ end and the in-read index sequence identifies their sample origin. deseq is used to estimate differential transcript abundance between conditions, independent of gene annotation. while the strand-specific tc 3′ ends can be linked to any gene annotation, here we use the ensembl gene build. fig. 2d shows the total number of regions called as peaks, the subset associated with a gene where deseq has estimated an adjusted p-value, genes showing differential transcript abundance between conditions with an adjusted p-value <= and finally the subset of the latter showing differential transcript abundance with a fold change > the analysis pipeline utilises a single configuration file describing the samples, the location of the sequence files and the conditions, and with one command converts duplicate-marked bam files into the detct pipeline output gene list tables. these tables show the closest strand-specific ensembl gene to the tc 3′ end, the region called as a peak, the unadjusted p-value, the adjusted p-value, the fold change between conditions and the count data. it is presented in tab or comma-separated tables or as an html table.fig. detct pipeline workflow. between nine and pairs of mutant and normal zebrafish embryos were collected from one clutch and rna extracted. a following dnasei treatment and chemical fragmentation, molecules representing the 3′ end of transcripts were enriched by pulldown using an anchored biotinylated oligo dt primer attached to streptavidin magnetic beads. an rna oligo matching part of the illumina read adapter was ligated onto the 5′ end, the rna eluted and annealed to an oligo comprising partial read illumina adapter followed by random bases, then an eight base indexing sequence specific to each sample and finally a base anchored polyt sequence. after reverse transcription the illumina adapter sequences were completed during library amplification. libraries were quantified, pooled in equimolar amounts and sequenced by illumina hiseq b after decoding the indexing sequence, the trimmed zebrafish sequences were mapped to the reference genome and duplicate reads were flagged. c the coordinate representing the transcript counting 3′ end was predicted using the base immediately 3′ of the polyt sequence in read. after calling peaks using all mapped read 2s the resulting counts were associated with their respective sample. the count data were used to identify differential transcript abundance between conditions using deseq and reported as a fold change with an adjusted p-value. the tc 3′ ends were matched to the closest ensembl transcript 3′ ends on the same strand. gene list tables were produced and ordered by the lowest adjusted p-value. these gene lists were filtered for genes showing differential transcript abundance using the adjusted p-value and the proximity of the tc 3′ end and ensembl gene end fig. library features. four sets of libraries were prepared from zebrafish embryo lines carrying a disrupted gene. morphologically abnormal and normal sibling embryos were collected from in-crossed lines and total rna extracted. for each mutant line to samples of indexed libraries were made and sequenced by illumina hiseq a the number of reads and the number mapping to the zv reference genome per library are shown. the total reads, the mean per library and the standard deviation are shown on the right. b for each library the proportion of reads identified as duplicates using outer mapping coordinates alone are shown on the x-axis and after accounting for the unique molecule identifier on the y-axis. c the reads were passed through the detct analysis pipeline. the number of read 2s mapped and the number of counts called under peaks as discrete regions are shown per library sample. the total reads, the mean per library and the standard deviation are shown on the right. d using the gene list output from the detct pipeline the chart shows the number of discrete regions identified per collection of libraries, the number of genes with an adjusted p-value from deseq, the number of genes with an adjusted p-value <= and the number of genes with an adjusted p-value <= plus a fold change  >  these analyses were performed with a stringent or relaxed proximity filter between the ensembl transcript and tc 3′ end data interpretation the experimental rationale suggests that all transcript counting 3′ ends should match the 3′ ends of transcripts. however, gene annotation is sometimes incomplete and occasionally both the annotation and tc 3′ end can arise from experimental artefact. therefore the detct output specifies the distance between the tc 3′ end and the nearest ensembl transcript 3′ end on the same strand. naturally, this assumes the transcript to be annotated correctly and this is not necessarily true. the ensembl transcript ends that exactly match the tc 3′ ends can be easily filtered from the results table. non-exact matches suggest incomplete gene annotation or novel alternative transcript ends and both situations can be validated by individual inspection. choosing a close proximity filter of the coordinates, such as between - bases and + bases, reduces the likelihood of a false positive match and using these criteria we were able to detect a mean of genes per experiment. in contrast a more relaxed proximity filter, such as between - to + identifies many more genes, but also finds more false positive ends. this is discussed further in the comparison to rna-seq below. one cause of false positive tc 3′ ends is oligo dt priming from polya or degenerate polya sequence within rna molecules or possibly from residual dna. we identified falsely primed tc 3′ ends during the detct analysis pipeline by examining the bases 3′ of the tc 3′ end and removed those which potentially derived from non-polya tail priming using the criteria described in the methods. one zebrafish-specific example is in the mitochondrial genome where a region rich in adenine in the rrna ensdarg accounts for between % to % of the total counts in the four test experiments. this example escaped our current polya filtering method. to mitigate false positive 3′ end calls, we have begun to prepare a list of tc 3′ ends we believe not to be true transcript tc 3′ ends. similarly, we are able to catalogue all the true positive tc 3′ ends we find and build cross-experiment profiles of regularly identified tc 3′ ends. we further filtered the results file by restricting the regions used to those with the 3′ coordinate of the region not more than bases upstream of the tc 3′ end. assessing the detct method to assess the variation between libraries we extracted rna from a pool of zebrafish embryos and made replicate transcript counting libraries with different indexing sequences using μg each. libraries were pooled, sequenced on an illumina miseq and analysed using the detct pipeline. the number of normalised counts for each genomic region called by detct for each library was determined. the regions were filtered for a maximum bp distance to an ensembl gene 3′ end, as described in additional file and compared using a pearson correlation. these libraries show our method displays good technical reproducibility.fig. technical replicate. twelve replicate transcript counting libraries were prepared from μg of a pool of wild-type zebrafish embryo total rna sample. libraries were sequenced by illumina miseq and analysed using the detct pipeline. a the normalised counts for each region were extracted. after filtering the data for only the regions which we would use to call differential transcript abundance the counts from all libraries were compared using a pearson correlation. cells coloured yellow in the pearson correlation are the most highly correlated while those in blue are the least correlated with a colour gradient inbetween. b in addition four concentrations of ercc spike mix were added in triplicate to the same libraries prior to library construction. we added the quantity suggested by the manufacturer, five times the quantity, one fifth of the quantity and one tenth of the quantity. the reads were mapped to the zebrafish reference sequence and ercc spike reference sequence. the diagram shows the ercc spikes represented by a circle in descending order of spike copy number in the mix on the x-axis and spike abundance on the y-axis. the blue circles show spikes identified in the detct pipeline while those in red were not found. c the detct analysis pipeline was run using six libraries at a time with three replicates as one pair of conditions and in all six possible condition combinations. the mean log fold change was calculated for all the spikes detected by the detct analysis and plotted against the expected log fold change as circles. each circle is labelled with the conditions being compared and observed log fold change over the expected log fold change. the numbers in brackets indicate how many spikes show differential transcript abundance to examine the method’s performance in quantifying differential transcript abundance we added ercc spike mix to the total rna prior to making the technical replicate libraries described above. we added the quantity of spike mix recommended by the manufacturer, five times the quantity, a fifth of the quantity and a tenth of the quantity, all in triplicate, to create four conditions. after analysis through the detct pipeline we identified spikes. the spikes with high copy number in the mix were detected, but not those with lower copy number. the mean log fold changes were calculated for the spikes detected in all combinations of condition and compared to the log fold change expected. we found a good correlation between observed and expected log fold change. additionally, in three of the pairwise comparisons no false positives were detected and in the other three only a total of six non-ercc regions were found to be differentially expressed, suggesting that the method has high specificity. our method was designed to maximise the number of tissue samples we can process with relatively shallow sequencing whilst still obtaining sufficient information to implicate gene networks modified by the condition change. the ability to make numerous transcript counting libraries has two main advantages. firstly, problem libraries resulting from sample loss or showing low complexity can be excluded from analysis and their removal has little effect on the statistical power of the analysis. similarly, samples with incorrectly assigned condition have less influence on the final result. secondly, increasing the number of replicates improves the statistical power of the analysis. to assess the impact of increasing the number of replicates whilst retaining the same total amount of sequencing we performed a permutation test using the zmp_ph data. reads were combined to produce two samples comprising, mutant and, wild-type reads. these were randomly split ten times into ten collections of pseudo-samples each containing equal numbers of reads but with the number of pairs of pseudo-samples ranging from two pairs to pairs. after passing these collections of pseudo-samples through the detct pipeline we identified transcript counting 3′ ends in all simulations. although this simulation has removed the variance from the original biological replicates it shows that increasing the number of libraries at the expense of read depth improves detection of tc 3′ ends showing differential abundance, as was previously noted.fig. sample number. the mutant and wild-type zmp_ph bam files were downsampled and merged to create two bam files, one from the mutant samples and one from the wild-type samples. the real index sequences in each bam file were then replaced with fake index sequences in order to simulate assigning the reads at random to between to different pseudo-samples. each of these permutations was repeated times, to give simulations in total, each of which was run through the detct pipeline. the chart shows the average total number of counts for each pseudo-sample and the number of tc 3′ ends with a relaxed ensembl transcript proximity filter showing differential transcript abundance with the mean above each group of ten many mrna expression pipelines use whole transcript rna-seq protocols and a range of analysis tools. we don’t intend to replace these methods but present an alternative method for high-throughput mrna expression screening. however, it is useful to compare the results of both methods. to this end we made detct libraries and standard non-directional polya pulldown illumina rna-seq libraries for two alleles. the same three wild-type and three mutant zebrafish total rna samples were processed for each method, plus two or six additional libraries using the detct protocol. we sequenced one hiseq lane equivalent for each allele by each method. potential duplicate reads were identified and eliminated. the umi in detct allowed a more accurate identification of duplicates and hence fewer reads were dismissed as duplicates. read 2s from the rna-seq data that mapped to the genome were compared to ensembl gene annotation to produce count data for each gene and the detct pipeline was used to extract count data linked to ensembl transcripts from the transcript counting reads. deseq was run on both sets of count data. even after removing duplicate reads there are generally more counts in the rna-seq count data. this is probably due to a drop in read quality following the oligo dt sequence in transcript counting read rna-seq initially identifies more genes, however, the gap between the methods is reduced substantially when regions with a low mean count are filtered out by deseq. the number of tc 3′ ends is higher than that of detct genes which suggests alternative 3′ ends. however, some represent false positive tc 3′ ends which escape our filter for false oligo dt priming. to assess the genes showing differential transcript abundance the full gene list was filtered for protein-coding genes with an adjusted p-value of <= and an absolute fold change >= between mutant and wild-type. rna-seq identifies more genes, but using a less stringent proximity filter between the tc 3′ end and the 3′ end of an ensembl transcript increases the detection rate in detct. for a direct comparison between the two methods we identified the genes with an adjusted p-value <= and absolute fold change >= for both methods and applied the less stringent proximity filter to the tc 3′ ends. the fold change for these genes in rna-seq and tc was compared and showed a good correlation with r2 = , which suggests the two methods are finding the same alterations in transcript abundance. we next looked at genes which show an adjusted p-value <= and an absolute fold change >= by one method, but failed to meet one or both criteria by the other. for both methods genes have an absolute fold change >= but fail to have sufficient power to call an adjusted p-value <= by one method. similarly, for genes where one method fails to show an absolute fold change >= the actual fold change is just below suggesting further genes where the two methods give comparable fold change results. we then applied the stringent tc 3′ end proximity filter to the same data which led to the removal of genes. examining the tc 3′ end of these genes showed they fell into two groups, either true tc 3′ ends or false tc 3′ ends assumed to be derived from experimental artefact. where both methods gave an adjusted p-value <= and an absolute fold change >= all were shown to be true ends. by contrast, in the gene sets only called by one method 11/ tc 3′ ends lost to more stringent filtering were false positives. together this analysis shows the removal of genes by increasing the stringency of the detct proximity filter resulted in losing true positives, but prevented calling false positive tc 3′ ends.table comparison of rna-seq and transcript counting lamc1 sa379 mdn1 sa1349 fig. comparing rna-seq and transcript counting. the genes identified as showing differential transcript abundance in experiment lamc1 sa with an adjusted p-value <= and an absolute fold change of >= by one or both methods at the less stringent proximity filter between the transcript counting 3′ end and the ensembl transcript end were plotted with their log fold change from rna-seq against the log fold change from tc. genes with an adjusted p-value <= and absolute fold change >= in both rna-seq and tc are shown as blue circles. genes which fail on one or both criteria in rna-seq but not in tc are shown as green circles and vice versa as red circles. a the log fold changes of genes with the tc 3′ end and ensembl transcript proximity filter between - and + the arrows highlight examples of genes which are not seen in graph b. b the log fold changes of genes with the proximity filter between - and + c a table showing the number of genes represented by each circle colour. the fourth and fifth columns show the genes which are lost during the more stringent proximity filtering and whether these are considered true positives or false positives after examining the tc 3′ end. note that for two genes where the rna-seq and tc data correlated, the closest tc 3′ end was found to be false, but a true end was identified further downstream and indicated by square brackets. the number in curly brackets indicates genes with a true tc 3′ end and an absolute fold change > but which fail to show an adjusted p-value <= in tc  we present a quantitative mrna transcript profiling package that starts with tissue samples and produces a gene list by counting the 3′ end of any polyadenylated transcripts using illumina sequencing. unlike whole transcript rna-seq each transcript is counted only once giving a more representative estimate of transcript abundance. short, rare transcripts are as likely to be represented in the sequence as long, rare transcripts. assaying only the 3′ end of transcripts is also more resilient to degraded rna samples, particularly if cells or tissues are compromised by a treatment. differential transcript abundance is identified based on genome sequence and is independent of gene annotation. this highlights regions containing novel transcripts as well as previously undescribed alternative transcript 3′ ends, which are implicated in biological processes related to the condition. each alternative transcript 3′ end is represented by discrete count data and has the potential to add layers of functional annotation to sequence at the 3′ ends of gene models. our streamlined library preparation reduces material loss and allows us to produce large numbers of libraries. working with more replicates not only provides more power to the differential expression calculation but also lessens the impact of occasional sample loss or failure. this allows us to measure transcriptional changes even in rare tissues within the whole organism. the addition of a unique molecular identifier helps us to assess the quality of each library preparation and we can remove any which underperform from the analysis. simultaneously processing a large number of samples allows us to ask complex questions, such as what is the impact of different concentrations of a compound at varying stages of embryonic development? similarly, we could screen large numbers of mutagenised individuals following an infection challenge to identify variation in infection response and inflammation. in both cases relatively shallow sequencing of each sample will result in reduced sensitivity, but can highlight a critical combination of conditions to explore further using deeper sequencing of selected libraries to increase the number of transcripts detected or by an alternative method such as rna-seq. we are also able to assay large numbers of wild-type samples producing a gene profile baseline that can be refined with data from each individual biological replicate. for precious or difficult to obtain rna samples conventional full transcript rna-seq methods are more suited, but for rapid assessment of many samples with a desire to implicate a molecular process, transcript counting is the better option. the gene lists provided by detct are the link between a living organism and the array of gene ontology, gene expression and gene interaction network data in the public domain. we believe the power of the differential expression transcript counting technique lies in effectively and efficiently bridging this gap. 
multiplex sequencing of transcript 3′ ends identifies differential transcript abundance independent of gene annotation. this method can be implemented on polyadenylated rna from any organism with an annotated reference genome and in any laboratory with access to illumina sequencing.
 we present a genome-wide messenger rna sequencing technique that converts small amounts of rna from many samples into molecular phenotypes. it encompasses all steps from sample preparation to sequence analysis and is applicable to baseline profiling or perturbation measurements.  multiplex sequencing of transcript 3′ ends identifies differential transcript abundance independent of gene annotation. we show that increasing biological replicate number while maintaining the total amount of sequencing identifies more differentially abundant transcripts.  this method can be implemented on polyadenylated rna from any organism with an annotated reference genome and in any laboratory with access to illumina sequencing. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords mrna transcript profilingrna-seqmolecular phenotypeissue-copyright-statement© the author 2015 