BACKGROUND
a variety of techniques, most notably tiling arrays  <cit> , massive sequencing of complimentary dnas   <cit>  and cdna tags  <cit> , has consistently identified the presence of rna molecules encoded from regions of the genome not currently annotated as exons of protein-coding rnas in human cells  <cit> .
collectively, these rna molecules were dubbed 'dark matter in the genome'  <cit> .
while these approaches underscored the complexity of transcriptional output from the human genome  <cit> , tiling arrays could predict only the complexity, but not the relative mass of 'dark matter' rna as a fraction of total transcription and, hence, they could not make strong conclusions about its importance.
similarly, the majority of methods for cdna sequencing have either used polya+ selected rna and/or amplification for library construction, both of which selectively omit significant amounts of rna and, hence, cannot accurately reflect the profile of all cellular rnas.
these technological limitations in assessing all cellular rnas have led to the notion that non-coding rnas  represent only a minor proportion of the total mass of rna molecules relative to protein-coding species and may arise only as a by-product of transcriptional noise  <cit> .
if so, they would not be worthy of additional analysis while others have suggested that the situation merits further study  <cit> .
the advent of next-generation sequencing technologies and their application to rna analysis  <cit>  provides the opportunity for us to re-assess the relative mass of the nc-rnas and, thus, directly pose the question of how prevalent and important this class of rnas might be.
however, care must be taken to ensure the rna analysed is representative of the entire cell and that the limitations of the technologies being used are understood to minimize the impact of any biases.
one initial report  <cit>  has suggested that nc-rnas represent only a minor fraction of rna in a mammalian cell and that most of these non-coding species consisted of intronic by-products generated from unspliced rnas.
in order to further address this question, we have used single-molecule sequencing  to reduce the potential for representational biases by avoiding amplification and minimizing sample preparation  <cit> .
we show that when total rna is sequenced rather than just the highly selected polyadenylated rnas, 'dark matter' rnas represent a significant fraction, sometimes even the majority, of the stable non-ribosomal, non-mitochondrial rnas  in a cell.
furthermore, gene-desert regions, currently devoid of annotations, are found to be expressed at abundant levels and rna species from these regions may be over-represented during or following neoplastic transformation.
we provide estimates of the relative mass of the 'dark matter' rna in cells of two organisms - human and fly.
to our knowledge, while used in a number of reports, this term has not been defined and, thus, we offer our definition of this term: 'dark matter' rna includes any rna whose function we do not currently understand and any rna that is currently un-annotated.
the latter category may include both protein-coding and nc-rnas.
a more detailed breakdown of what is included and how these species might arise is provided in the additional file 1: definition of dark matter.
another parameter requiring definition is the measurement for transcriptome complexity used here.
in general, the fraction of 'dark matter' rna could be measured either as a fraction of nucleotide bases in the genome that are used to encode these transcripts or as the relative mass of these rnas as a fraction of mass of all rnas in the cell after removing ribosomal and mtrnas.
most techniques, notably tiling arrays have provided the former estimate.
importantly, however, this estimator alone cannot provide a complete description of the complexity of 'dark matter' rna as it is possible that it can be highly complex in terms of the sequences present in this population and yet comprise a minute fraction of total rna population by mass.
this has been suggested to be the case by the recent report of van bakel et al.
who estimate that 'dark matter' rna constitutes 12% of all rna by mass  <cit> .
in that report, the authors have leveraged the ability of a next-generation sequencing to count individual sequencing reads and assign them to either known or 'dark matter' rna categories.
the fraction of total reads in each category of rna represents its relative mass.
we have also used the relative mass as defined by fraction of reads falling within a certain category as an estimator of the abundance of rna.
we used the number of reads as opposed to the total number of nucleotides covered by the reads to estimate the relative mass of rna since the sms reads do not have the same lengths.
RESULTS
in order to evaluate the true nature of the complete cellular rna population, we analysed rna from various tissue sources and different rna preparations in two different species: human and drosophila .
the filtered sequence reads used in this work have been deposited in the national center for biotechnology information short read archive .
polya+ rna, total rna and total rna depleted of ribosomal rna were studied.
as described in the materials and methods section, rna preparations were treated with saturating concentrations of dnase i followed by purification over columns to eliminate small amounts of any remaining partially-digested small dna pieces.
this step can also remove a substantial fraction of rnas below  <dig> bases and, thus, this work is primarily aimed at rnas longer than  <dig> bases.
the priming of rna with random hexamers and reverse transcription generated first-strand cdna which was then tailed with terminal transferase and deoxyadenosine triphosphate .
this was followed by sequencing using single-molecule sequencing after hybridization to a flow cell surface with covalently bound poly dt <dig>  <cit> .
no amplification, ligation or size selection were used in cdna preparation, minimizing methodological biases  <cit> .
each cdna sample was sequenced on one or more channels of a heliscope genetic analysis system with the resulting filtered reads  aligned to the complete human or drosophila genomes supplemented with the sequence of the complete ribosomal dna  repeat .
reverse transcription is known to produce spurious second-strand cdna  <cit> , so reads from both strands of the genome were combined in order to avoid sense/antisense artifacts.
the unique mapping reads from human tissue sources were further filtered to exclude sequences aligned to rdna sequences, the mitochondrial genome, as well as to genomic repeats annotated by the repeatmasker program as rrna.
after filtering, the remaining informative reads were used for subsequent analyses, including a comparison to the known annotations defined by the university of california, santa cruz  genes or the flybase genes tracks from the ucsc browser.
as shown in figure  <dig> and additional file 2: table s <dig>  the highest fraction of non-exonic reads consistently occurred in total rna  and total rna depleted of ribosomal rna .
strikingly lower proportions of non-exonic reads  occurred in the polya+ fraction, although they still accounted for approximately 19%-33% of the informative reads.
two different polya+ rna preparation methods for human liver rna, using either magnetic beads or oligo-dt cellulose purification, produced similar yields of informative reads  that overlap known exons.
tissue source had a modest effect on the fraction of non-exonic rnas, with the normal brain having the highest proportion of ncrna transcripts, followed by the leukaemia cell line k <dig> and then the liver, which showed the lowest level.
since these data were not biased by potential effects of amplification or ligation, they are more likely to accurately represent the relative mass of different populations of rnas than other techniques which require additional processing steps.
any skewing in favour of non-exonic rnas that could be induced by the process of rna purification from the cell or reverse transcription would have also occurred to a similar extent in other studies.
overall, this suggests that the un-annotated transcriptome, as defined by the reads mapping outside exons of protein-coding rnas, represents the majority of rnas by mass, comprising up to two-thirds of all non-ribosomal, non-mtrna in a human cell.
however, of interest and consistent with the previous observation based on tiling arrays  <cit> , rna from adult drosophila flies did not have as much non-exonic rnas by mass as human cells, comprising approximately 11% of the polya+ rna and 19% of the ribominus rna.
this is consistent with earlier observations based on tiling arrays suggesting that drosophila does not have as much 'dark matter' rna as human cells  <cit>  and supports the previous suggestion that the fraction of the 'dark matter' rna increases with increasing organismal complexity  <cit> .
as with the human samples, the fly samples have a higher fraction of 'dark matter' rna in total rather than polya+ rna .
the data and analysis below refer to human samples unless otherwise specified.
in order to further investigate the complexity of transcriptional output in human samples, we have analysed six different ewing family of tumours  samples from four different patients.
these included a pair of matching cell lines derived from a single patient, one from the primary-tumour established in culture prior to chemotherapy  and one from a metastatic-lesion cultured after chemotherapy .
a second matched set of samples was derived from primary and metastatic eft tumours from a different patient.
an additional primary and metastatic tumour sample were each obtained from different patients .
the primary and metastatic tumour samples served to exclude the possibility that the observed non-exonic rnas are an artifact of tumour cell line culture conditions.
figure  <dig> and additional file 3: table s <dig> summarizes the analysis of these samples profiled using ribominus rna.
consistent with the previous results, non-exonic rnas represented 43%-63% of all non-ribosomal, non-mtrnas by mass.
primary tumour samples had more non-exonic rnas  than the cell lines , which could be due to the more heterogeneous nature of the tumour tissues.
tumours are known to be comprised of multiple cell types with normal tissue stroma yielding up to 50% of the total rna mass and different cell-types could have different non-exonic rnas expression profiles  <cit> .
thus, overall, 30%-52% of all informative reads in the human ribominus rna did not overlap exons and fell within intronic regions  with an average of 40%, representing a significantly higher fraction than the  <dig> % reported by van bakel et al.
<cit> .
intronic reads could simply represent unspliced rnas and to some extent this will be the case.
however, the abundance and specificity that we find among intronic rnas indicate they are not simply bystanders but suggest a more complex role.
introns have been shown in several instances to harbour known non-coding functional rnas, such as the kcnq1ot <dig> transcript involved in imprinting  <cit> .
by rna mass in a human cell, transcripts emanating from intronic sequences approximately equal that of exonic sequences but this large amount of intronic sequence cannot be explained just by the fact that introns are longer and, thus, accumulate more reads.
the density of reads from individual introns can be quite abundant and similar to, or higher than, that of exonic regions.
this is exemplified by the known ncrna kcnq1ot <dig> embedded within the protein-coding kcnq <dig> locus and transcribed from the opposite strand, indicating it is not simply a splicing artifact .
additional examples in loci not currently known to harbour ncrnas are shown on figure 4b.
in contrast, rna transcripts from some introns are virtually absent .
in the immp2l and slco1b3/lst-3tm <dig> loci , the density of reads in intronic regions is equal to, or higher than, in the exonic regions, with a notable enrichment of intronic sequences in the ribominus compared to polya+ rna.
overall, such examples are not rare.
if we consider all spliced ucsc genes annotations expressed in k <dig> defined by the presence of at least  <dig> reads per  <dig> m in intronic and exonic regions, we identify  <dig>  transcripts.
of these,  <dig> have at least one intron whose average read density is equal to, or higher than, the density of the exons.
it is also important to emphasize that we obtained these results from an analysis of cellular total rna, not nuclear-enriched rnas and, thus, the relative mass of introns in the total rna pool is impressive.
overall, our results are consistent with the fact that introns produce abundant stand-alone transcripts.
consistent with this, louro et al.
<cit>  have found that over 80% of human refseq loci have  <dig>  transcriptional units formed by clusters of expressed sequence tags  that map totally within introns, suggesting that introns contain transcripts that are not part of pre-messenger rna.
importantly, these units are not randomly distributed in the genome  <cit> .
furthermore, intronic rnas were previously found to represent a compilation of different rna species in addition to independent sense and antisense transcripts - transcripts retaining specific introns as well as rna molecules with alternative exons and alternative un-translated regions  <cit> .
although less frequent than intronic reads, a substantial number of ribominus rna reads also mapped to intergenic space.
intergenic reads comprise 24%-37% of all non-exonic reads and approximately 15%-18% of all informative reads in the ribominus rna samples .
interestingly, we also noticed a number of very long, abundant intergenic transcribed regions spanning hundreds of kbs of genome virtually devoid of annotations, as illustrated in figures 4c and 4d.
similar to the intronic rnas, these transcripts were also enriched in the ribominus rnas, although still detectable in polya+ rna samples.
these very long intergenic rnas were found at significantly greater abundance in cell lines derived from malignancies, such as k <dig>  compared to rna from the normal liver and brain .
in order to further investigate this in cancer, we analysed the six eft samples and found many very long intergenic transcribed regions specific to these eft samples, one example of which appears in figure 5a.
a region of approximately  <dig> kb on chromosome  <dig> was observed in total rna from the eft samples, but not in rnas from k <dig>  normal liver or the brain.
interestingly, transcription from this region occurs in at least three of the six eft samples, including a cell line and two solid tumour samples, suggesting that it is not an artifact of in vitro cell culture.
furthermore, we detected long, intergenic transcribed regions specific to k <dig> cells that were transcribed from an approximately  <dig> kb intergenic region on chromosome  <dig> and not found in eft samples .
while any given genomic region could be rearranged, amplified or deleted in the cancerous cells, it is unlikely that the same regions are rearranged in samples derived from different sources as exemplified by the eft samples derived from four different patients.
in order to identify all such expressed intergenic regions, we pooled non-genic reads from the eft samples and k <dig> cells and defined intergenic rna domains of at least  <dig> kb in length .
overall,  <dig> such domains were found .
we thus propose to classify these intergenic regions as very long intergenic non-coding  regions.
as described below, expression from vlinc regions can be substantial and significantly higher in cancerous tissues than in normal tissues.
the total stable rna output of a vlinc region can be quite substantial based on relative mass.
in the case of k <dig>  if both the vlinc and known protein-coding transcripts are ranked by relative mass, the top expressed vlinc  has a higher relative mass than the ninth most expressed known transcript which corresponds to the eukaryotic translation elongation factor  <dig>  eef <dig> .
of further interest is the observation that seven of the top eight expressed annotated transcripts in this cell line correspond to annotated nc-rnas, such as malat <dig> and rmrp .
thus, the only protein-coding rna that is higher in relative mass than the top vlinc rnas is an isoform of epsilon globin gene .
the total stable rna output of vlinc regions tends to be higher in cancerous tissues than in normal ones.
the maximum number of normalized counts for the seven tumour  and two normal samples  was calculated for each vlinc .
the tumour/normal fold ratio of the maximal counts was then calculated for each vlinc .
the median ratio was approximately  <dig> , showing that most of the vlincs are more abundant in the tumour samples .
furthermore, of the  <dig> vlinc regions,  <dig> had a ratio that was double or more and, for  <dig> vlincs, this ratio was  <dig> times higher.
in contrast, there were three and zero vlinc regions whose maximal read count in a normal tissue was two and  <dig> times higher, respectively, than that in a tumour tissue/cell line.
additionally, the total stable rna output of vlinc regions can differentiate between different tumours and even tumours of the same category derived from different patients.
a non-parametric spearman correlation analysis among all of the pair-wise combinations of the eft , normal brain and liver sample and the k <dig> leukaemia cell line was performed based on the expression values for vlinc and ucsc genes .
the following average correlation metrics based on all relevant pair-wise correlations were then calculated: between the eft samples from the same patient; between the eft samples from different patients; between the eft and the normal tissues; and between the eft samples and the leukaemia cell line.
the correlation matrix is shown in additional file 6: table s <dig>  for the ucsc genes, sum of expression within each annotated transcript of the ucsc genes database was calculated across all samples .
the transcripts were then sorted by the sum expression level and the top two quartiles  were used to calculate the correlations.
this was done to avoid inflation of the correlations due to transcripts not expressed in any of the tissues.
the average correlation between the eft samples from the same patient was, as expected, quite high using either the expression of ucsc genes or vlincs .
it is worth noting that the pair of eft samples from the same patient were represented by the primary and metastatic tumour isolated at significantly different time periods  from each other and, in the case of the one pair represented by the chla-9/- <dig> cell lines, cultured independently.
this suggests that the expression of vlinc regions can be maintained for significant periods of time and suggests an important biological role.
the average correlation between the eft samples from the different patients was still high using expression of ucsc genes but lower using vlincs:  <dig>  and  <dig> , respectively.
the latter was similar to those calculated between eft samples and normal tissues:  <dig>  and  <dig>  based on ucsc genes and vlinc regions, respectively.
interestingly, expression levels of vlincs in eft tumours were extremely different from the leukaemia k <dig> cell line .
however, the average correlation between the two tumour types using known genes was still relatively high:  <dig> .
although the number of tumour tissues used in this study is low, these observations suggest that the 'dark matter' rnas produced from the vlinc regions can differentiate between different tumours, perhaps even better than the annotated genes.
the ability to differentiate molecularly defined subgroups in the case of the eft and, probably other tumours, is especially interesting since tumours are known to be heterogeneous and, in the case of eft, are often referred to as a 'family of tumours'.
further, the diagnosis of eft is often challenging due to its non-variant, microscopic appearance  <cit> .
this data further confirms the notion of heterogeneity and suggests that the 'dark matter' rnas could, in fact, be as good as or even better a discriminator than other methods, including coding rnas, which are not unique to any given tumour, for identifying molecular subgroups within this family of tumours or to distinguish it from other tumour types.
previously, human large intergenic non-coding  rnas were identified in normal human embryonic and stem cell lines  <cit>  and we sought to determine whether the link regions were overlapping.
in fact, the majority of the link transcribed regions we have identified did not overlap the known human lacuna regions and, thus, represent novel rnas that are also large, interagency and non-coding, as exemplified in the four examples shown .
these latter regions have known lacuna regions located nearby, without overlap, while the former do not have lacuna regions in their vicinity.
furthermore, the interagency regions identified here achieve much greater lengths than known lacunas, with a median size of ~ <dig> kb versus  <dig> kb for the lacunas .
overall, 37%  of the regions overlapped the k4-k <dig> domains harbouring lacunas as reported by khalid et al.
<cit> .
however, even when overlapping, the lacuna regions corresponded to only a fraction of our interagency regions: the percent of base pairs in the interagency regions found here only overlapped the lacuna regions by approximately 19% .
however, the overlap between the two categories of the interagency transcribed regions is highly significant .
this is interesting because the two regions were detected using different experimental approaches and even more so, they were found in very different types of tissues.
links were originally found in normal embryonic tissues/cell lines, while links were originally found in tumour tissues.
this suggests a tantalizing possibility that links are also expressed very early in development, perhaps at very specific stages of development, which explains why they so far have eluded detection and the determination of function during early development.
we also analysed whether vlincs could be conserved across species, at least in terms of the syntenic location.
in order to do this we took advantage of the data denoting 5' and 3' ends of mouse transcripts  from the fantom <dig> consortium  <cit> .
we took the coordinates of  <dig>  such clusters and converted these clusters into regions of the hg <dig> version of the human genome .
overall,  <dig>  clusters could be remapped to hg <dig>  of the  <dig>  clusters,  <dig> were intergenic and, of those,  <dig> spanned at least  <dig> kb of genomic space - the minimal size of a vlinc region.
these  <dig> clusters are thus analogous to our vlinc regions based on size and location in the intergenic regions.
overall,  <dig> vlinc regions were found to overlap  <dig>  fantom <dig> clusters by at least  <dig> kb .
sometimes a vlinc region overlapped more than one cluster .
overall, given the small span of the genomic space covered by the vlincs and the  <dig> fantom <dig> clusters, the overlap was significant .
furthermore, in the cases of  <dig> vlincs, the region of overlap was more the 50% of the length of either the vlinc or the fantom <dig> cluster it overlapped .
thus, some vlincs are conserved in mouse at least in the syntenic locations.
the fact that only  <dig>  vlincs are conserved in mouse is probably due to the fact that:  only  <dig> such regions were identified in fantom <dig>  probably due to usage of polya+ rna for library constructions by the fantom <dig> consortium, significantly lowering potential for overlap; and  differences in biological material - vlincs appear to be highly tissue specific.
it is difficult to currently estimate the complexity and properties of rna molecules harboured within the vlinc regions.
what is currently annotated explains a minute fraction of the complexity found in these regions, both in terms of the base pair coverage and fraction of reads, even if we include additional databases, such as refseq, ensemble and all spliced ests.
for example, exons of the spliced ests cover only ~ <dig> % of the genomic sequences covered by all vlincs:  <dig>  out  <dig> , <dig> bp correspondingly.
also, taking k <dig> as an example, approximately 2% of all reads  that fall within the vlinc regions can be accounted for by the exons of the spliced ests.
this number climbs to  <dig> % if exons of all spliced ests, ucsc genes, refseq and ensemble are combined .
as mentioned above, however, when summarized over each vlinc region, these transcripts are quite abundant based on how many reads one needs to minimally cover all non-ribosomal, non-mt rnas in a human cell.
nevertheless, elucidation of the structures of transcripts within the vlinc regions is of paramount importance.
however, this can only be done by targeting each region and cloning and sequencing each full-length cdna, to the extent that this is even possible.
both cloning and 5' and 3' race  experiments may prove impossible if the individual transcripts turn out to be very long, necessitating new strategies for characterizing these regions.
given that we sequenced total as well as polya-selected rna, we can also estimate the relative mass of non-ribosomal, non-mt rnas in both fractions and also estimate the number of reads required to fully sequence this fraction from one cell.
sequencing of total rna not depleted of ribosomal rna yields a value of about 6%  of non-ribosomal and non-mt rna as a fraction of total rna in cell by mass.
based on the amount of rna present in one mammalian cell  and an average read of  <dig> bases , we would need greater than  <dig> million reads to sequence all non-ribosomal, non-mt rna in a single human cell.
since the number of reads in this study of ribominus rnas represents only a fraction of the total number of reads estimated above, we can safely assume that we sequenced less than one cellular equivalent of rna.
thus, the readily apparent detection of the above non-exonic rnas, as shown in figures 3- <dig>  suggests that their expression occurs at significant levels compared to protein encoding rnas and certainly much higher than one copy per cell.
we can also estimate that the polya-selected material probably represents only a small fraction  of non-ribosomal, non-mtrna.
this estimate follows from the typical yield of a single polya selection, which represents approximately 1%-2% of starting total rna amounts, of which approximately 50% is still ribosomal or mtrna and, from the observation made in this study, that total rna contains on the order of 6%  non-ribosomal and non-mt rna.
CONCLUSIONS
in summary, we show that the relative mass of 'dark-matter' rnas is both substantial and, in some cases, greater than that of protein encoding transcripts.
we find that the 'dark-matter' rna fraction typically represents 50%-65% of all non-ribosomal, non-mt rna by mass in many normal and neoplastic tissues.
we derived these estimates using single-molecule sequencing of cdnas that avoids amplification, ligation and cdna size-selection steps, all of which are known to increase the potential for bias prior to sequencing  <cit> .
this number is probably an underestimate since we considered only reads that fall within non-exonic regions.
exonic sequences that are part of transcripts that have no known function were not counted due to the limitations discussed earlier.
also, because the issues with reverse-transcriptase mentioned above make it difficult to unambiguously define the strand of all rnas, antisense exonic reads would be counted as exonic reads, also decreasing the fraction of the true 'dark matter' rnas.
the minimal sample preparation employed in this study more closely reflects the native rna population abundance compared with amplification-based, next-generation platforms.
we also find that total cellular rna has a much higher complexity than polya+ rna, which is highly enriched for transcripts comprised of currently annotated exonic sequences across a variety of cell types and tissues.
a methodological choice that assumes rnas need to be polyadenylated in order to be worthy of study by discarding everything without a polya tail is, thus, certain to leave gaping holes in understanding of the transcriptome.
we further show that the intergenic space in regions devoid of known, annotated transcripts harbours very long stretches , of abundantly transcribed and likely non-protein encoding rna species, even though we cannot exclude a possibility that some of transcripts in the vlinc regions may encode proteins.
our understanding of the repertoire of these types of rnas remains far from complete, as we have identified hundreds of new loci by profiling just seven samples from two types of tumours.
furthermore, the statistically significant overlap between the intergenic vlinc rnas found in tumours in this work and the lincrnas found in the normal human embryonal and stem cells  <cit>  raises the interesting hypothesis that these rnas in the vlinc regions function in specific capacities early in development, only to be silenced in differentiated cells and reactivated during malignant transformation.
these transcripts appear to be readily identifiable in tumours, particularly poorly differentiated, 'stem cell like' tumours such as eft, but less so in normal tissues, which suggests functional roles in both development and oncogenesis.
such expressed rna species may also have potential usefulness as diagnostic biomarkers for specific cancer types, as exemplified previously by the hotair lincrna  <cit>  and other non-coding rnas  <cit>  as well as in other diseases  <cit> .
it would be interesting to determine whether vlincs share common functional properties with macro ncrna identified as regulators of imprinting  <cit> .
further analysis of more tumours from a variety of individuals will be required in order to establish conclusive connections between expressed regions and cancer aetiology.
materials and methods
rna
total rna was obtained from commercial sources for the k <dig> cell line  and normal liver and normal brain tissues .
total rna was isolated from a pair of cell lines from a patient with eft, one established at diagnosis  and the second from a nodal metastasis obtained after chemotherapy   <cit> .
most eft tumours express a fusion protein  generated from a translocation and fusion between ewsr <dig> and an ets family gene.
the fusion protein acts as a transcription regulator and transcriptome sequencing of this tumour is probably to reveal mis-expressed rnas involved in neoplastic transformation.
cell lines tested mycoplasma-free.
cell line identity was validated using patient bone marrow by the short tandem repeat  assay  <cit>  and verified by str at time of nucleic acid extraction for these experiments.
cell lines can be obtained from the children's oncology group repository at http://www.cogcell.org.
the eft tissues were obtained from the repository at the children's hospital of los angeles, usa, and the detail of the specimens could be obtained from dr timothy triche .
briefly, their properties were as follows:
<dig>  primary no.
<dig> tumour: sacral tumour; age - <dig> years and  <dig> months; ewsr <dig> gene rearrangement positive by fish
<dig>  metastasis no.
<dig> tumour: lung metastasis; same patient as primary no.
<dig> tumour
<dig>  primary no.
<dig> tumour: paraspinal tumour; age -  <dig> years and  <dig> months; ewsr <dig> gene rearrangement positive by fish
<dig>  metastasis no.
<dig> tumour: metastasis of a primary thoraco-abdominal tumour; age - <dig> years and  <dig> months; real-time polymerase chain reaction and sequence analysis confirmed type i fusion of ewsr1/fli- <dig> genes.
before further fractionation, total rna was treated with dnase i as follows:  <dig> μg of total rna was mixed with  <dig> μl of 10x dnase i buffer ;  <dig> μl of rnaseout ; and  <dig> μl of recombinant dnase i  and incubated for  <dig> min at 37°c.
the rna was then purified using the rneasy minelute kit .
the dnase i-treated total rna was either unfractionated  or fractionated using one of the following methods:  depleted of ribosomal rna  using the ribominus kit ;  polya+ fraction was selected using a magnetic bead-based purification kit ; or,  polya+ fraction was selected using the oligo-dt cellulose method purist kit, ambion, texas, usa).
preparation of rna for sequencing
except where noted, 100- <dig> ng of dnase i -treated rna was mixed with the following reagents from the superscript iii kit .
first,  <dig> μl of  <dig> ng/μl random hexamers and  <dig> μl of  <dig> mm dntps were added in the total volume of  <dig> μl.
the mixture was then placed in a thermocycler and heat denatured at 65°c for  <dig> min followed by rapid cooling on ice.
next,  <dig> μl of 10x cdna synthesis buffer,  <dig> μl of  <dig>  m dtt and  <dig> μl of  <dig> mm mgcl <dig> were added.
the samples were returned to the thermocycler and allowed to incubate at 15°c for  <dig> min.
then,  <dig>  μl of rnaseout and  <dig>  μl of superscript iii reverse transcriptase were added and the samples incubated at 25°c for  <dig> min, 42°c for  <dig> min, 55°c for  <dig> min and 70°c for  <dig> min.
after reverse transcription, rna was removed by adding  <dig> μl of rnaseh  and  <dig> μl of rnase if  to each sample and incubating at 37°c for  <dig> min.
the cdna was then purified by two rounds of purification over performa columns  and quantified using a nanodrop spectrophotometer.
next, a 3' poly-a tail was added to the cdna samples.
cdna  was mixed with a control oligo to monitor tail length and water in a total volume of  <dig>  μl.
the mixture was denatured at 95°c for  <dig> min followed by rapid cooling on ice.
five microlitres of  <dig>  mm cocl <dig>   <dig> μl of 10x terminal deoxynucleotidyl transferase  buffer ,  <dig> μl of  <dig> μm datp and  <dig>  μl of tdt  was then added and the samples incubated at 42°c for  <dig> h and at 70°c for  <dig> min.
the 3' ends of the poly-a-tailed cdna were then blocked with biotin-ddatp.
the sample was denatured at 95°c for  <dig> min followed by rapid cooling on ice.
we then added  <dig>  μl of  <dig> mm biotin-ddatp  and  <dig>  μl of tdt followed by incubation at 37°c for  <dig> min and 70°c for  <dig> min.
the control oligo was removed by digestion with the user enzyme .
then  <dig> μl of the user enzyme  was added to the sample and incubated at 37°c for  <dig> min.
the sample was then purified using ampure beads  by bringing the volume up to  <dig> μl with water and adding  <dig> μl of the ampure beads followed by incubation at room temperature for  <dig> min with agitation.
the beads were then captured on a magnetic stand and washed twice with 70% ethanol.
the beads were allowed to air dry for 5- <dig> min, resuspended in  <dig> μl of water and left open for  <dig> min on the magnet.
the eluate was collected and the beads were resuspended again in  <dig> μl of water and then left for  <dig> min on the magnet.
the eluate was collected again and combined with the first eluate.
typically, samples were hybridized to the heliscope flow cell in  <dig> μl at a loading concentration of 100- <dig> pm.
sms reads to the genome
sms reads were trimmed for leading t homopolymers and were filtered for reads with a minimal length of  <dig> bases after trimming using a suite of helicos tools available at http://open.helicosbio.com/mwiki/index.php/releases and described at http://open.helicosbio.com/helisphere_user_guide/index.html.
alignments were conducted with indexdpgenomic software freely available on the helicos website http://open.helicosbio.com/mwiki/index.php/releases.
the aligner maximizes the aligned yield of sms reads due to the ability to align reads in which the predominant error is represented by deletions common in single-molecule sequencing.
for the genomic alignments, reads were aligned to the ncbiv <dig> version of the genome supplemented with the complete ribosomal repeat unit  or the dm <dig> version of the fly genome supplemented with the complete ribosomal repeat unit  using the helicos basic analysis pipeline http://open.helicosbio.com/helisphere_user_guide/index.html.
the sequence reads were filtered to include reads with a minimal length of  <dig> bases and aligned using a stringent normalized score of  <dig>   and  <dig>  .
aligned reads were further filtered for reads having a unique best alignment score.
the mapped reads were then compared to ucsc genes and flybase gene annotations to partition them into exonic, intronic and intergenic reads.
the normalized score was defined as follows:
score = /length_read
for example, in the following alignment:
tag sequence cctccgtgttgttccagcc-cagtgctcgcagg
ref sequence c-tccgtgttgttccagccacagtgctcgcagg
length of alignment block: 33
length of tag sequence: 32
number of matches: 31
number of errors: 2
score:  -  =  <dig> -  <dig> = 147
normalized score = 147/ <dig> =  <dig>
identification of the vlinc rna domains
reads not overlapping the genic regions from the ewing tumours and k <dig> were pooled and used to generate the densities of read coverage throughout the genome.
the vlinc domains were defined in the integrated genome browser http://www.bioviz.org/igb/download.shtml by applying the following thresholds to the density graphs for each chromosome:
threshold = 80th percentile
maxgap = 5000
minrun = 50000
comparison with the fantom  <dig> dataset
the coordinates of the  <dig>  reliable 5'/3' clusters described in the table  <dig> of the fantom <dig> consortium paper  <cit>  were downloaded from the fantom <dig> website at ftp://fantom.gsc.riken.jp/fantom3/boundary_set/pair53_clusters.txt.gz.
the coordinates were lifted from the mm <dig> version of the mouse genome to hg <dig> in a two step process using the ucsc liftover tool: mm5- > mm <dig> and mm8- > hg <dig>
abbreviations
eft: ewing family of tumours; fish: fluorescent in situ hybridization; lincnas: large intergenic non-coding; mtrna: mitochondrial rna; ncrna: non-coding rna; rdna: recomninant dna; vlinc: very long intergenic non-coding; sms: single molecule sequencing; str: short tandem repeat; tdt: terminal deoxynucleotidyl transferase.
authors' contributions
pk, pm, jft and tjt designed the study.
gsl, cpr and tjt provided the experimental material.
pk performed the analysis.
pk, gsl, rja, jft and tjt wrote the manuscript and provided input on the analysis and the presentation of the data.
pk, tr and fo performed the molecular biological experiments.
phbs and gr provided intellectual contribution to the manuscript.
all authors read and approved the final manuscript.
supplementary material
additional file 1
definition of the 'dark matter' rna.
provides a definition of what kind of rna molecules should be included in the 'dark matter' rna realm and argumentation for their inclusion.
click here for file
additional file 2
table s <dig>  distribution of single-molecule sequencing reads obtained from polya+, rrna-depleted  and total rna from different tissues among different types of annotations.
click here for file
additional file 3
table s <dig>  distribution of single-molecule sequencing reads obtained from ribosomal rna-depleted  rna from ewing family of tumours solid tumours and cell lines.
click here for file
additional file 4
table s <dig>  coordinates of the very long intergenic regions and their normalized expression levels in each tissue.
click here for file
additional file 5
table s <dig>  normalized expression levels of university of california santa cruz genes in each tissue.
the total number of reads within each of the regions was normalized to  <dig> m non-ribosomal non-mitochondrial unique reads in each sample.
click here for file
additional file 6
table s <dig>  spearman correlation matrix based on expression of very long intergenic non-codings and university of california santa cruz genes between different tumour and normal tissues
click here for file
additional file 7
table s <dig>  overlap between very long intergenic non-codings and fantom <dig> intergenic clusters of over  <dig> kb.
click here for file
acknowledgements
tjt wishes to acknowledge and thank dr jim jacobson, acting associate director and chief, diagnostic biomarker and technologies branch, cancer diagnosis program, dctd nci  for his continuous support of the work that led to this study as an outgrowth of nci funded uo <dig> studies .
we thank dr jason farrar for his careful reading of the manuscript; drs sylvain foissac, doron lipson and stan letovsky for their helpful discussions; drs michele wing and yuxia zhan in the genome core at chla for their help with the rna preparation; dr hiro shimada for help with the eft tissue handling; and members of the triche laboratory for their help.
the work was supported by a grant from the department of defense , w81xwh-07-1- <dig> and through support of the children's oncology group  chair's grant, u- <dig> ca <dig> from the nih; by rett's initiative and by the uo <dig> specs program .
rja is partially supported through the endowed king fahd professorship in pediatric oncology at johns hopkins.
this publication was made possible in part by grant rc <dig> hg <dig>  from the national human genetics research institute  at the national institutes of health.
its contents are solely the responsibility of the authors and do not necessarily represent the official views of nhgri.

 a variety of techniques, most notably tiling arrays, massive sequencing of complimentary dnas and cdna tags, has consistently identified the presence of rna molecules encoded from regions of the genome not currently annotated as exons of protein-coding rnas in human cells. collectively, these rna molecules were dubbed'dark matter in the genome'. while these approaches underscored the complexity of transcriptional output from the human genome, tiling arrays could predict only the complexity, but not the relative mass of'dark matter' rna as a fraction of total transcription and, hence, they could not make strong conclusions about its importance. similarly, the majority of methods for cdna sequencing have either used polya+ selected rna and/or amplification for library construction, both of which selectively omit significant amounts of rna and, hence, cannot accurately reflect the profile of all cellular rnas. these technological limitations in assessing all cellular rnas have led to the notion that non-coding rnas represent only a minor proportion of the total mass of rna molecules relative to protein-coding species and may arise only as a by-product of transcriptional noise. if so, they would not be worthy of additional analysis while others have suggested that the situation merits further study. the advent of next-generation sequencing technologies and their application to rna analysis provides the opportunity for us to re-assess the relative mass of the nc-rnas and, thus, directly pose the question of how prevalent and important this class of rnas might be. however, care must be taken to ensure the rna analysed is representative of the entire cell and that the limitations of the technologies being used are understood to minimize the impact of any biases. one initial report has suggested that nc-rnas represent only a minor fraction of rna in a mammalian cell and that most of these non-coding species consisted of intronic by-products generated from unspliced rnas. in order to further address this question, we have used single-molecule sequencing to reduce the potential for representational biases by avoiding amplification and minimizing sample preparation. we show that when total rna is sequenced rather than just the highly selected polyadenylated rnas,'dark matter' rnas represent a significant fraction, sometimes even the majority, of the stable non-ribosomal, non-mitochondrial rnas in a cell. furthermore, gene-desert regions, currently devoid of annotations, are found to be expressed at abundant levels and rna species from these regions may be over-represented during or following neoplastic transformation. we provide estimates of the relative mass of the'dark matter' rna in cells of two organisms - human and fly. to our knowledge, while used in a number of reports, this term has not been defined and, thus, we offer our definition of this term:'dark matter' rna includes any rna whose function we do not currently understand and any rna that is currently un-annotated. the latter category may include both protein-coding and nc-rnas. a more detailed breakdown of what is included and how these species might arise is provided in the additional file 1: definition of dark matter. another parameter requiring definition is the measurement for transcriptome complexity used here. in general, the fraction of'dark matter' rna could be measured either as a fraction of nucleotide bases in the genome that are used to encode these transcripts or as the relative mass of these rnas as a fraction of mass of all rnas in the cell after removing ribosomal and mtrnas. most techniques, notably tiling arrays have provided the former estimate. importantly, however, this estimator alone cannot provide a complete description of the complexity of'dark matter' rna as it is possible that it can be highly complex in terms of the sequences present in this population and yet comprise a minute fraction of total rna population by mass. this has been suggested to be the case by the recent report of van bakel et al. who estimate that'dark matter' rna constitutes 12% of all rna by mass. in that report, the authors have leveraged the ability of a next-generation sequencing to count individual sequencing reads and assign them to either known or'dark matter' rna categories. the fraction of total reads in each category of rna represents its relative mass. we have also used the relative mass as defined by fraction of reads falling within a certain category as an estimator of the abundance of rna. we used the number of reads as opposed to the total number of nucleotides covered by the reads to estimate the relative mass of rna since the sms reads do not have the same lengths.  in order to evaluate the true nature of the complete cellular rna population, we analysed rna from various tissue sources and different rna preparations in two different species: human and drosophila. the filtered sequence reads used in this work have been deposited in the national center for biotechnology information short read archive. polya+ rna, total rna and total rna depleted of ribosomal rna were studied. as described in the materials and methods section, rna preparations were treated with saturating concentrations of dnase i followed by purification over columns to eliminate small amounts of any remaining partially-digested small dna pieces. this step can also remove a substantial fraction of rnas below bases and, thus, this work is primarily aimed at rnas longer than bases. the priming of rna with random hexamers and reverse transcription generated first-strand cdna which was then tailed with terminal transferase and deoxyadenosine triphosphate. this was followed by sequencing using single-molecule sequencing after hybridization to a flow cell surface with covalently bound poly dt. no amplification, ligation or size selection were used in cdna preparation, minimizing methodological biases. each cdna sample was sequenced on one or more channels of a heliscope genetic analysis system with the resulting filtered reads aligned to the complete human or drosophila genomes supplemented with the sequence of the complete ribosomal dna repeat. reverse transcription is known to produce spurious second-strand cdna, so reads from both strands of the genome were combined in order to avoid sense/antisense artifacts. the unique mapping reads from human tissue sources were further filtered to exclude sequences aligned to rdna sequences, the mitochondrial genome, as well as to genomic repeats annotated by the repeatmasker program as rrna. after filtering, the remaining informative reads were used for subsequent analyses, including a comparison to the known annotations defined by the university of california, santa cruz genes or the flybase genes tracks from the ucsc browser. as shown in figure and additional file 2: table s the highest fraction of non-exonic reads consistently occurred in total rna and total rna depleted of ribosomal rna. strikingly lower proportions of non-exonic reads occurred in the polya+ fraction, although they still accounted for approximately 19%-33% of the informative reads. two different polya+ rna preparation methods for human liver rna, using either magnetic beads or oligo-dt cellulose purification, produced similar yields of informative reads that overlap known exons. tissue source had a modest effect on the fraction of non-exonic rnas, with the normal brain having the highest proportion of ncrna transcripts, followed by the leukaemia cell line k and then the liver, which showed the lowest level. since these data were not biased by potential effects of amplification or ligation, they are more likely to accurately represent the relative mass of different populations of rnas than other techniques which require additional processing steps. any skewing in favour of non-exonic rnas that could be induced by the process of rna purification from the cell or reverse transcription would have also occurred to a similar extent in other studies. overall, this suggests that the un-annotated transcriptome, as defined by the reads mapping outside exons of protein-coding rnas, represents the majority of rnas by mass, comprising up to two-thirds of all non-ribosomal, non-mtrna in a human cell. however, of interest and consistent with the previous observation based on tiling arrays, rna from adult drosophila flies did not have as much non-exonic rnas by mass as human cells, comprising approximately 11% of the polya+ rna and 19% of the ribominus rna. this is consistent with earlier observations based on tiling arrays suggesting that drosophila does not have as much'dark matter' rna as human cells and supports the previous suggestion that the fraction of the'dark matter' rna increases with increasing organismal complexity. as with the human samples, the fly samples have a higher fraction of'dark matter' rna in total rather than polya+ rna. the data and analysis below refer to human samples unless otherwise specified. in order to further investigate the complexity of transcriptional output in human samples, we have analysed six different ewing family of tumours samples from four different patients. these included a pair of matching cell lines derived from a single patient, one from the primary-tumour established in culture prior to chemotherapy and one from a metastatic-lesion cultured after chemotherapy. a second matched set of samples was derived from primary and metastatic eft tumours from a different patient. an additional primary and metastatic tumour sample were each obtained from different patients. the primary and metastatic tumour samples served to exclude the possibility that the observed non-exonic rnas are an artifact of tumour cell line culture conditions. figure and additional file 3: table s summarizes the analysis of these samples profiled using ribominus rna. consistent with the previous results, non-exonic rnas represented 43%-63% of all non-ribosomal, non-mtrnas by mass. primary tumour samples had more non-exonic rnas than the cell lines, which could be due to the more heterogeneous nature of the tumour tissues. tumours are known to be comprised of multiple cell types with normal tissue stroma yielding up to 50% of the total rna mass and different cell-types could have different non-exonic rnas expression profiles. thus, overall, 30%-52% of all informative reads in the human ribominus rna did not overlap exons and fell within intronic regions with an average of 40%, representing a significantly higher fraction than the % reported by van bakel et al.. intronic reads could simply represent unspliced rnas and to some extent this will be the case. however, the abundance and specificity that we find among intronic rnas indicate they are not simply bystanders but suggest a more complex role. introns have been shown in several instances to harbour known non-coding functional rnas, such as the kcnq1ot transcript involved in imprinting. by rna mass in a human cell, transcripts emanating from intronic sequences approximately equal that of exonic sequences but this large amount of intronic sequence cannot be explained just by the fact that introns are longer and, thus, accumulate more reads. the density of reads from individual introns can be quite abundant and similar to, or higher than, that of exonic regions. this is exemplified by the known ncrna kcnq1ot embedded within the protein-coding kcnq locus and transcribed from the opposite strand, indicating it is not simply a splicing artifact. additional examples in loci not currently known to harbour ncrnas are shown on figure 4b. in contrast, rna transcripts from some introns are virtually absent. in the immp2l and slco1b3/lst-3tm loci, the density of reads in intronic regions is equal to, or higher than, in the exonic regions, with a notable enrichment of intronic sequences in the ribominus compared to polya+ rna. overall, such examples are not rare. if we consider all spliced ucsc genes annotations expressed in k defined by the presence of at least reads per m in intronic and exonic regions, we identify transcripts. of these, have at least one intron whose average read density is equal to, or higher than, the density of the exons. it is also important to emphasize that we obtained these results from an analysis of cellular total rna, not nuclear-enriched rnas and, thus, the relative mass of introns in the total rna pool is impressive. overall, our results are consistent with the fact that introns produce abundant stand-alone transcripts. consistent with this, louro et al. have found that over 80% of human refseq loci have transcriptional units formed by clusters of expressed sequence tags that map totally within introns, suggesting that introns contain transcripts that are not part of pre-messenger rna. importantly, these units are not randomly distributed in the genome. furthermore, intronic rnas were previously found to represent a compilation of different rna species in addition to independent sense and antisense transcripts - transcripts retaining specific introns as well as rna molecules with alternative exons and alternative un-translated regions. although less frequent than intronic reads, a substantial number of ribominus rna reads also mapped to intergenic space. intergenic reads comprise 24%-37% of all non-exonic reads and approximately 15%-18% of all informative reads in the ribominus rna samples. interestingly, we also noticed a number of very long, abundant intergenic transcribed regions spanning hundreds of kbs of genome virtually devoid of annotations, as illustrated in figures 4c and 4d. similar to the intronic rnas, these transcripts were also enriched in the ribominus rnas, although still detectable in polya+ rna samples. these very long intergenic rnas were found at significantly greater abundance in cell lines derived from malignancies, such as k compared to rna from the normal liver and brain. in order to further investigate this in cancer, we analysed the six eft samples and found many very long intergenic transcribed regions specific to these eft samples, one example of which appears in figure 5a. a region of approximately kb on chromosome was observed in total rna from the eft samples, but not in rnas from k normal liver or the brain. interestingly, transcription from this region occurs in at least three of the six eft samples, including a cell line and two solid tumour samples, suggesting that it is not an artifact of in vitro cell culture. furthermore, we detected long, intergenic transcribed regions specific to k cells that were transcribed from an approximately kb intergenic region on chromosome and not found in eft samples. while any given genomic region could be rearranged, amplified or deleted in the cancerous cells, it is unlikely that the same regions are rearranged in samples derived from different sources as exemplified by the eft samples derived from four different patients. in order to identify all such expressed intergenic regions, we pooled non-genic reads from the eft samples and k cells and defined intergenic rna domains of at least kb in length. overall, such domains were found. we thus propose to classify these intergenic regions as very long intergenic non-coding regions. as described below, expression from vlinc regions can be substantial and significantly higher in cancerous tissues than in normal tissues. the total stable rna output of a vlinc region can be quite substantial based on relative mass. in the case of k if both the vlinc and known protein-coding transcripts are ranked by relative mass, the top expressed vlinc has a higher relative mass than the ninth most expressed known transcript which corresponds to the eukaryotic translation elongation factor eef. of further interest is the observation that seven of the top eight expressed annotated transcripts in this cell line correspond to annotated nc-rnas, such as malat and rmrp. thus, the only protein-coding rna that is higher in relative mass than the top vlinc rnas is an isoform of epsilon globin gene. the total stable rna output of vlinc regions tends to be higher in cancerous tissues than in normal ones. the maximum number of normalized counts for the seven tumour and two normal samples was calculated for each vlinc. the tumour/normal fold ratio of the maximal counts was then calculated for each vlinc. the median ratio was approximately, showing that most of the vlincs are more abundant in the tumour samples. furthermore, of the vlinc regions, had a ratio that was double or more and, for vlincs, this ratio was times higher. in contrast, there were three and zero vlinc regions whose maximal read count in a normal tissue was two and times higher, respectively, than that in a tumour tissue/cell line. additionally, the total stable rna output of vlinc regions can differentiate between different tumours and even tumours of the same category derived from different patients. a non-parametric spearman correlation analysis among all of the pair-wise combinations of the eft, normal brain and liver sample and the k leukaemia cell line was performed based on the expression values for vlinc and ucsc genes. the following average correlation metrics based on all relevant pair-wise correlations were then calculated: between the eft samples from the same patient; between the eft samples from different patients; between the eft and the normal tissues; and between the eft samples and the leukaemia cell line. the correlation matrix is shown in additional file 6: table s for the ucsc genes, sum of expression within each annotated transcript of the ucsc genes database was calculated across all samples. the transcripts were then sorted by the sum expression level and the top two quartiles were used to calculate the correlations. this was done to avoid inflation of the correlations due to transcripts not expressed in any of the tissues. the average correlation between the eft samples from the same patient was, as expected, quite high using either the expression of ucsc genes or vlincs. it is worth noting that the pair of eft samples from the same patient were represented by the primary and metastatic tumour isolated at significantly different time periods from each other and, in the case of the one pair represented by the chla-9/- cell lines, cultured independently. this suggests that the expression of vlinc regions can be maintained for significant periods of time and suggests an important biological role. the average correlation between the eft samples from the different patients was still high using expression of ucsc genes but lower using vlincs: and, respectively. the latter was similar to those calculated between eft samples and normal tissues: and based on ucsc genes and vlinc regions, respectively. interestingly, expression levels of vlincs in eft tumours were extremely different from the leukaemia k cell line. however, the average correlation between the two tumour types using known genes was still relatively high:. although the number of tumour tissues used in this study is low, these observations suggest that the'dark matter' rnas produced from the vlinc regions can differentiate between different tumours, perhaps even better than the annotated genes. the ability to differentiate molecularly defined subgroups in the case of the eft and, probably other tumours, is especially interesting since tumours are known to be heterogeneous and, in the case of eft, are often referred to as a'family of tumours'. further, the diagnosis of eft is often challenging due to its non-variant, microscopic appearance. this data further confirms the notion of heterogeneity and suggests that the'dark matter' rnas could, in fact, be as good as or even better a discriminator than other methods, including coding rnas, which are not unique to any given tumour, for identifying molecular subgroups within this family of tumours or to distinguish it from other tumour types. previously, human large intergenic non-coding rnas were identified in normal human embryonic and stem cell lines and we sought to determine whether the link regions were overlapping. in fact, the majority of the link transcribed regions we have identified did not overlap the known human lacuna regions and, thus, represent novel rnas that are also large, interagency and non-coding, as exemplified in the four examples shown. these latter regions have known lacuna regions located nearby, without overlap, while the former do not have lacuna regions in their vicinity. furthermore, the interagency regions identified here achieve much greater lengths than known lacunas, with a median size of ~ kb versus kb for the lacunas. overall, 37% of the regions overlapped the k4-k domains harbouring lacunas as reported by khalid et al.. however, even when overlapping, the lacuna regions corresponded to only a fraction of our interagency regions: the percent of base pairs in the interagency regions found here only overlapped the lacuna regions by approximately 19%. however, the overlap between the two categories of the interagency transcribed regions is highly significant. this is interesting because the two regions were detected using different experimental approaches and even more so, they were found in very different types of tissues. links were originally found in normal embryonic tissues/cell lines, while links were originally found in tumour tissues. this suggests a tantalizing possibility that links are also expressed very early in development, perhaps at very specific stages of development, which explains why they so far have eluded detection and the determination of function during early development. we also analysed whether vlincs could be conserved across species, at least in terms of the syntenic location. in order to do this we took advantage of the data denoting 5' and 3' ends of mouse transcripts from the fantom consortium. we took the coordinates of such clusters and converted these clusters into regions of the hg version of the human genome. overall, clusters could be remapped to hg of the clusters, were intergenic and, of those, spanned at least kb of genomic space - the minimal size of a vlinc region. these clusters are thus analogous to our vlinc regions based on size and location in the intergenic regions. overall, vlinc regions were found to overlap fantom clusters by at least kb. sometimes a vlinc region overlapped more than one cluster. overall, given the small span of the genomic space covered by the vlincs and the fantom clusters, the overlap was significant. furthermore, in the cases of vlincs, the region of overlap was more the 50% of the length of either the vlinc or the fantom cluster it overlapped. thus, some vlincs are conserved in mouse at least in the syntenic locations. the fact that only vlincs are conserved in mouse is probably due to the fact that: only such regions were identified in fantom probably due to usage of polya+ rna for library constructions by the fantom consortium, significantly lowering potential for overlap; and differences in biological material - vlincs appear to be highly tissue specific. it is difficult to currently estimate the complexity and properties of rna molecules harboured within the vlinc regions. what is currently annotated explains a minute fraction of the complexity found in these regions, both in terms of the base pair coverage and fraction of reads, even if we include additional databases, such as refseq, ensemble and all spliced ests. for example, exons of the spliced ests cover only ~ % of the genomic sequences covered by all vlincs: out, bp correspondingly. also, taking k as an example, approximately 2% of all reads that fall within the vlinc regions can be accounted for by the exons of the spliced ests. this number climbs to % if exons of all spliced ests, ucsc genes, refseq and ensemble are combined. as mentioned above, however, when summarized over each vlinc region, these transcripts are quite abundant based on how many reads one needs to minimally cover all non-ribosomal, non-mt rnas in a human cell. nevertheless, elucidation of the structures of transcripts within the vlinc regions is of paramount importance. however, this can only be done by targeting each region and cloning and sequencing each full-length cdna, to the extent that this is even possible. both cloning and 5' and 3' race experiments may prove impossible if the individual transcripts turn out to be very long, necessitating new strategies for characterizing these regions. given that we sequenced total as well as polya-selected rna, we can also estimate the relative mass of non-ribosomal, non-mt rnas in both fractions and also estimate the number of reads required to fully sequence this fraction from one cell. sequencing of total rna not depleted of ribosomal rna yields a value of about 6% of non-ribosomal and non-mt rna as a fraction of total rna in cell by mass. based on the amount of rna present in one mammalian cell and an average read of bases, we would need greater than million reads to sequence all non-ribosomal, non-mt rna in a single human cell. since the number of reads in this study of ribominus rnas represents only a fraction of the total number of reads estimated above, we can safely assume that we sequenced less than one cellular equivalent of rna. thus, the readily apparent detection of the above non-exonic rnas, as shown in figures 3- suggests that their expression occurs at significant levels compared to protein encoding rnas and certainly much higher than one copy per cell. we can also estimate that the polya-selected material probably represents only a small fraction of non-ribosomal, non-mtrna. this estimate follows from the typical yield of a single polya selection, which represents approximately 1%-2% of starting total rna amounts, of which approximately 50% is still ribosomal or mtrna and, from the observation made in this study, that total rna contains on the order of 6% non-ribosomal and non-mt rna.  in summary, we show that the relative mass of'dark-matter' rnas is both substantial and, in some cases, greater than that of protein encoding transcripts. we find that the'dark-matter' rna fraction typically represents 50%-65% of all non-ribosomal, non-mt rna by mass in many normal and neoplastic tissues. we derived these estimates using single-molecule sequencing of cdnas that avoids amplification, ligation and cdna size-selection steps, all of which are known to increase the potential for bias prior to sequencing. this number is probably an underestimate since we considered only reads that fall within non-exonic regions. exonic sequences that are part of transcripts that have no known function were not counted due to the limitations discussed earlier. also, because the issues with reverse-transcriptase mentioned above make it difficult to unambiguously define the strand of all rnas, antisense exonic reads would be counted as exonic reads, also decreasing the fraction of the true'dark matter' rnas. the minimal sample preparation employed in this study more closely reflects the native rna population abundance compared with amplification-based, next-generation platforms. we also find that total cellular rna has a much higher complexity than polya+ rna, which is highly enriched for transcripts comprised of currently annotated exonic sequences across a variety of cell types and tissues. a methodological choice that assumes rnas need to be polyadenylated in order to be worthy of study by discarding everything without a polya tail is, thus, certain to leave gaping holes in understanding of the transcriptome. we further show that the intergenic space in regions devoid of known, annotated transcripts harbours very long stretches, of abundantly transcribed and likely non-protein encoding rna species, even though we cannot exclude a possibility that some of transcripts in the vlinc regions may encode proteins. our understanding of the repertoire of these types of rnas remains far from complete, as we have identified hundreds of new loci by profiling just seven samples from two types of tumours. furthermore, the statistically significant overlap between the intergenic vlinc rnas found in tumours in this work and the lincrnas found in the normal human embryonal and stem cells raises the interesting hypothesis that these rnas in the vlinc regions function in specific capacities early in development, only to be silenced in differentiated cells and reactivated during malignant transformation. these transcripts appear to be readily identifiable in tumours, particularly poorly differentiated,'stem cell like' tumours such as eft, but less so in normal tissues, which suggests functional roles in both development and oncogenesis. such expressed rna species may also have potential usefulness as diagnostic biomarkers for specific cancer types, as exemplified previously by the hotair lincrna and other non-coding rnas as well as in other diseases. it would be interesting to determine whether vlincs share common functional properties with macro ncrna identified as regulators of imprinting. further analysis of more tumours from a variety of individuals will be required in order to establish conclusive connections between expressed regions and cancer aetiology. materials and methods rna total rna was obtained from commercial sources for the k cell line and normal liver and normal brain tissues. total rna was isolated from a pair of cell lines from a patient with eft, one established at diagnosis and the second from a nodal metastasis obtained after chemotherapy. most eft tumours express a fusion protein generated from a translocation and fusion between ewsr and an ets family gene. the fusion protein acts as a transcription regulator and transcriptome sequencing of this tumour is probably to reveal mis-expressed rnas involved in neoplastic transformation. cell lines tested mycoplasma-free. cell line identity was validated using patient bone marrow by the short tandem repeat assay and verified by str at time of nucleic acid extraction for these experiments. cell lines can be obtained from the children's oncology group repository at http://www.cogcell.org. the eft tissues were obtained from the repository at the children's hospital of los angeles, usa, and the detail of the specimens could be obtained from dr timothy triche. briefly, their properties were as follows: primary no. tumour: sacral tumour; age - years and months; ewsr gene rearrangement positive by fish metastasis no. tumour: lung metastasis; same patient as primary no. tumour primary no. tumour: paraspinal tumour; age - years and months; ewsr gene rearrangement positive by fish metastasis no. tumour: metastasis of a primary thoraco-abdominal tumour; age - years and months; real-time polymerase chain reaction and sequence analysis confirmed type i fusion of ewsr1/fli- genes. before further fractionation, total rna was treated with dnase i as follows: μg of total rna was mixed with μl of 10x dnase i buffer; μl of rnaseout; and μl of recombinant dnase i and incubated for min at 37°c. the rna was then purified using the rneasy minelute kit. the dnase i-treated total rna was either unfractionated or fractionated using one of the following methods: depleted of ribosomal rna using the ribominus kit; polya+ fraction was selected using a magnetic bead-based purification kit; or, polya+ fraction was selected using the oligo-dt cellulose method purist kit, ambion, texas, usa). preparation of rna for sequencing except where noted, 100- ng of dnase i -treated rna was mixed with the following reagents from the superscript iii kit. first, μl of ng/μl random hexamers and μl of mm dntps were added in the total volume of μl. the mixture was then placed in a thermocycler and heat denatured at 65°c for min followed by rapid cooling on ice. next, μl of 10x cdna synthesis buffer, μl of m dtt and μl of mm mgcl were added. the samples were returned to the thermocycler and allowed to incubate at 15°c for min. then, μl of rnaseout and μl of superscript iii reverse transcriptase were added and the samples incubated at 25°c for min, 42°c for min, 55°c for min and 70°c for min. after reverse transcription, rna was removed by adding μl of rnaseh and μl of rnase if to each sample and incubating at 37°c for min. the cdna was then purified by two rounds of purification over performa columns and quantified using a nanodrop spectrophotometer. next, a 3' poly-a tail was added to the cdna samples. cdna was mixed with a control oligo to monitor tail length and water in a total volume of μl. the mixture was denatured at 95°c for min followed by rapid cooling on ice. five microlitres of mm cocl μl of 10x terminal deoxynucleotidyl transferase buffer, μl of μm datp and μl of tdt was then added and the samples incubated at 42°c for h and at 70°c for min. the 3' ends of the poly-a-tailed cdna were then blocked with biotin-ddatp. the sample was denatured at 95°c for min followed by rapid cooling on ice. we then added μl of mm biotin-ddatp and μl of tdt followed by incubation at 37°c for min and 70°c for min. the control oligo was removed by digestion with the user enzyme. then μl of the user enzyme was added to the sample and incubated at 37°c for min. the sample was then purified using ampure beads by bringing the volume up to μl with water and adding μl of the ampure beads followed by incubation at room temperature for min with agitation. the beads were then captured on a magnetic stand and washed twice with 70% ethanol. the beads were allowed to air dry for 5- min, resuspended in μl of water and left open for min on the magnet. the eluate was collected and the beads were resuspended again in μl of water and then left for min on the magnet. the eluate was collected again and combined with the first eluate. typically, samples were hybridized to the heliscope flow cell in μl at a loading concentration of 100- pm. sms reads to the genome sms reads were trimmed for leading t homopolymers and were filtered for reads with a minimal length of bases after trimming using a suite of helicos tools available at http://open.helicosbio.com/mwiki/index.php/releases and described at http://open.helicosbio.com/helisphere_user_guide/index.html. alignments were conducted with indexdpgenomic software freely available on the helicos website http://open.helicosbio.com/mwiki/index.php/releases. the aligner maximizes the aligned yield of sms reads due to the ability to align reads in which the predominant error is represented by deletions common in single-molecule sequencing. for the genomic alignments, reads were aligned to the ncbiv version of the genome supplemented with the complete ribosomal repeat unit or the dm version of the fly genome supplemented with the complete ribosomal repeat unit using the helicos basic analysis pipeline http://open.helicosbio.com/helisphere_user_guide/index.html. the sequence reads were filtered to include reads with a minimal length of bases and aligned using a stringent normalized score of and. aligned reads were further filtered for reads having a unique best alignment score. the mapped reads were then compared to ucsc genes and flybase gene annotations to partition them into exonic, intronic and intergenic reads. the normalized score was defined as follows: score = /length_read for example, in the following alignment: tag sequence cctccgtgttgttccagcc-cagtgctcgcagg ref sequence c-tccgtgttgttccagccacagtgctcgcagg length of alignment block: 33 length of tag sequence: 32 number of matches: 31 number of errors: 2 score: - = - = 147 normalized score = 147/ = identification of the vlinc rna domains reads not overlapping the genic regions from the ewing tumours and k were pooled and used to generate the densities of read coverage throughout the genome. the vlinc domains were defined in the integrated genome browser http://www.bioviz.org/igb/download.shtml by applying the following thresholds to the density graphs for each chromosome: threshold = 80th percentile maxgap = 5000 minrun = 50000 comparison with the fantom dataset the coordinates of the reliable 5'/3' clusters described in the table of the fantom consortium paper were downloaded from the fantom website at ftp://fantom.gsc.riken.jp/fantom3/boundary_set/pair53_clusters.txt.gz. the coordinates were lifted from the mm version of the mouse genome to hg in a two step process using the ucsc liftover tool: mm5- > mm and mm8- > hg abbreviations eft: ewing family of tumours; fish: fluorescent in situ hybridization; lincnas: large intergenic non-coding; mtrna: mitochondrial rna; ncrna: non-coding rna; rdna: recomninant dna; vlinc: very long intergenic non-coding; sms: single molecule sequencing; str: short tandem repeat; tdt: terminal deoxynucleotidyl transferase. authors' contributions pk, pm, jft and tjt designed the study. gsl, cpr and tjt provided the experimental material. pk performed the analysis. pk, gsl, rja, jft and tjt wrote the manuscript and provided input on the analysis and the presentation of the data. pk, tr and fo performed the molecular biological experiments. phbs and gr provided intellectual contribution to the manuscript. all authors read and approved the final manuscript. supplementary material additional file 1 definition of the'dark matter' rna. provides a definition of what kind of rna molecules should be included in the'dark matter' rna realm and argumentation for their inclusion. click here for file additional file 2 table s distribution of single-molecule sequencing reads obtained from polya+, rrna-depleted and total rna from different tissues among different types of annotations. click here for file additional file 3 table s distribution of single-molecule sequencing reads obtained from ribosomal rna-depleted rna from ewing family of tumours solid tumours and cell lines. click here for file additional file 4 table s coordinates of the very long intergenic regions and their normalized expression levels in each tissue. click here for file additional file 5 table s normalized expression levels of university of california santa cruz genes in each tissue. the total number of reads within each of the regions was normalized to m non-ribosomal non-mitochondrial unique reads in each sample. click here for file additional file 6 table s spearman correlation matrix based on expression of very long intergenic non-codings and university of california santa cruz genes between different tumour and normal tissues click here for file additional file 7 table s overlap between very long intergenic non-codings and fantom intergenic clusters of over kb. click here for file acknowledgements tjt wishes to acknowledge and thank dr jim jacobson, acting associate director and chief, diagnostic biomarker and technologies branch, cancer diagnosis program, dctd nci for his continuous support of the work that led to this study as an outgrowth of nci funded uo studies. we thank dr jason farrar for his careful reading of the manuscript; drs sylvain foissac, doron lipson and stan letovsky for their helpful discussions; drs michele wing and yuxia zhan in the genome core at chla for their help with the rna preparation; dr hiro shimada for help with the eft tissue handling; and members of the triche laboratory for their help. the work was supported by a grant from the department of defense, w81xwh-07-1- and through support of the children's oncology group chair's grant, u- ca from the nih; by rett's initiative and by the uo specs program. rja is partially supported through the endowed king fahd professorship in pediatric oncology at johns hopkins. this publication was made possible in part by grant rc hg from the national human genetics research institute at the national institutes of health. its contents are solely the responsibility of the authors and do not necessarily represent the official views of nhgri.
we show that the relative mass of rna whose function and/or structure we do not understand, as a proportion of all non-ribosomal, non-mitochondrial human rna, can be greater than that of protein-encoding transcripts. discovery that the transcriptional output of the human genome is far more complex than predicted by the current set of protein-coding annotations and that most rnas produced do not appear to encode proteins has transformed our understanding of genome complexity and suggests new paradigms of genome regulation. we conclude that'dark matter' rna can constitute the majority of non-ribosomal, non-mitochondrial-rna and a significant fraction arises from numerous very long, intergenic transcribed regions that could be involved in neoplastic transformation. we further show the presence of a large number of very long, abundantly-transcribed regions in intergenic space and further show that expression of these regions is associated with neoplastic transformation. these overlap some regions found previously in normal human embryonic tissues and raises an interesting hypothesis as to the function of these ncrnas in both early development and neoplastic transformation.
 discovery that the transcriptional output of the human genome is far more complex than predicted by the current set of protein-coding annotations and that most rnas produced do not appear to encode proteins has transformed our understanding of genome complexity and suggests new paradigms of genome regulation. however, the fraction of all cellular rna whose function we do not understand and the fraction of the genome that is utilized to produce that rna remain controversial. this is not simply a bookkeeping issue because the degree to which this un-annotated transcription is present has important implications with respect to its biologic function and to the general architecture of genome regulation. for example, efforts to elucidate how non-coding rnas regulate genome function will be compromised if that class of rnas is dismissed as simply'transcriptional noise'.  we show that the relative mass of rna whose function and/or structure we do not understand, as a proportion of all non-ribosomal, non-mitochondrial human rna, can be greater than that of protein-encoding transcripts. this observation is obscured in studies that focus only on polya-selected rna, a method that enriches for protein coding rnas and at the same time discards the vast majority of rna prior to analysis. we further show the presence of a large number of very long, abundantly-transcribed regions in intergenic space and further show that expression of these regions is associated with neoplastic transformation. these overlap some regions found previously in normal human embryonic tissues and raises an interesting hypothesis as to the function of these ncrnas in both early development and neoplastic transformation.  we conclude that'dark matter' rna can constitute the majority of non-ribosomal, non-mitochondrial-rna and a significant fraction arises from numerous very long, intergenic transcribed regions that could be involved in neoplastic transformation. 