BACKGROUND
rna next generation sequencing  technology has replaced microarrays as the primary platform for gene expression profiling .
this same technology has also been used to study the expression of micrornas , although there are fewer direct comparisons between mirna arrays and mirna-seq .
during mirna-seq library preparation, small rnas are selected by electrophoresis on a size typically ranging from  <dig> to  <dig> nucleotides .
this range of size selection allows for the capture of many other species of small rnas in addition to mirnas .
one of the major byproducts of size selection for mirnas is trna-derived small rnas   <cit> .
this may be the result of either active cleavage or an artifact of the mirna-seq library construction.
full length or parent trnas are adaptor molecules composed of rna with a length typically ranging from  <dig> to  <dig> nts.
they serve as the physical link that translates messenger rna into protein.
cleavage of trnas by the rnase iii enzyme, angiogenin, may occur in a number of reactive conditions to produce trna-derived halves   <cit> .
likewise, trnas can also be cleaved in a dicer-dependent manner or as an in vitro phenomenon by incubation with mgcl <dig> or nuclease s <dig>  <cit> .
all of these processes result in cleavage within a hairpin loop.
d loop cleavage results in a 19 nt trna-derived fragment , and cleavage of the anticodon loop before the anti-codon affords a 33 nt trh  <cit> .
these 5’ tdr products, in particular the 33-base trh, may be captured by mirna-seq.
in addition, apparent tdrs and trnas with base errors found in rna-seq datasets may arise from the extensive chemical modifications of parent trna nucleotides.
while these modifications contribute physiologically to rna stability, function, and structure  <cit> , they may be either misincorporated during the requisite reverse transcription step  <cit> , resulting in specific base errors  <cit> , or result in blockage of reverse transcription, with production of a truncated cdna product  <cit> .
therefore, sequencing datasets may contain many reads that appear to be tdrs, but in fact are actually derived from full-length parent trna templates.
the quantification of trnas and tdrs warrant study through high-throughput sequencing.
aberrations in trnas have been linked to a range of diseases, including neurological disorders , cancer , type ii diabetes , and mitochondrial diseases .
in addition, tdrs have been identified in a number of cancer cell lines, and some tdrs may play defined biological roles, such as the role of trf- <dig> in proliferation  <cit> .
thus far, only one study has found a link between somatic mutations in mitochondrial trna , aka mt-tl1) and myelodysplastic syndromes   <cit> .
the potential functionality of trna in mds remains largely unexplored.
as part of our efforts to study mds, we have previously examined mirna expression in mds  <cit> .
mds is a common disease of the elderly   <cit>  characterized by ineffective maturation of hematopoietic cells.
this ineffective maturation manifests clinically as low blood counts and morphologically as dysplasia.
the mds state has been associated with expression changes in mirnas.
several entities have been discovered on 5q, which contribute to the disease phenotype in the 5q minus syndrome subtype of mds  <cit> .
mir- <dig> and mir- <dig> have recently been identified in the commonly deleted region, and their knockdown appears to be associated with the thrombocytosis and dysmegakaryopoiesis seen in 5q minus syndrome  <cit> .
researchers, including those in our group, have identified several other mirnas whose expression is dysregulated with the diagnosis of other subtypes of mds .
published research on other types of small rnas in mds is limited.
mds has a median survival of only 18– <dig> months  <cit> , with death resulting from either cytopenic complications or transformation to acute myeloid leukemia   <cit> .
currently, there are only two fda-approved drugs for the treatment of all types of mds, both of which are dna methyltransferase inhibitors : 5-azacytidine, and decitabine.
however, only 40–47 % of patients achieve the clinically meaningful responses of hematologic improvement or better with dnmtis  <cit> .
although several recent studies have suggested the use of mutational analysis by dna next generation sequencing   <cit> , genetics  <cit> , or the presence or absence of key pharmacodynamic and pharmacokinetic markers  <cit>  to predict response, as yet, there are no clinical or laboratory parameters in routine clinical practice that accurately predict response to dnmtis .
the methylation status of specific loci has been shown to be prognostically relevant in the treatment of a closely related disease, chronic myelomonocytic leukemia, by decitabine  <cit> .
in this study, we examined the expression of trnas in paired pre- and post-treatment samples from patients with the diagnosis of mds on receiving treatment with dnmtis.
in mirna sequencing data, the abundance of trna fragments vastly outnumbers mirnas.
this gives us a great opportunity to study trna expression using mirna sequencing data through their anticodon sequences.
methods
ethics statement
this study was conducted in accordance with the declaration of helsinki with the approval of the institutional review boards at roswell park cancer institute and vanderbilt university medical center.
consent statement
written informed consent for the patient-derived specimens was obtained prior to the acquisition in all cases.
the consents were approved by the institutional review boards at roswell park cancer institute and vanderbilt university medical center.
specimen description
bone marrow aspirate specimens   were obtained from the roswell park cancer institute with appropriate approval by the institutional review boards at roswell park and vanderbilt university medical center.
twenty-two of the  <dig> bone marrow samples were from control patients .
the other  <dig> samples were from mds patients.
out of  <dig> mds samples,  <dig> were pre-treatment specimens, and  <dig> were obtained after treatment with 5-azacytidine or decitabine, both dnmtis.
for each of the pre- and post-treatment paired samples, a response score  was obtained based on clinical and pathological criteria to indicate the effectiveness of treatment   <cit> .
bone marrow mononuclear cells  were isolated from fresh, unsorted bone marrow aspirate specimens  by the ficoll method using cellgro lymphocyte separation medium .
after performing a cell count, the cells were re-pelleted from dulbecco’s phosphate buffered saline solution and resuspended in gibco recovery freezing medium , and then frozen gradually, prior to storage in liquid nitrogen.
cells were frozen at a density of  <dig> to  <dig> million cells per ml.
rna isolation and rna sequencing
total rna  was isolated using a mirvana mirna isolation kit  per manufacturer’s instructions.
sequencing library construction was performed on the totrna from all  <dig> samples, each obtained from a single bone marrow aspirate specimen.
rnas were isolated by the mirvana rna isolation kit, which was found to perform better than the trizol mirna isolation kit  <cit> .
libraries were prepared using the truseq small rna sample preparation kit .
the small rna protocol specifically ligates rna adapters to mature mirnas that have a 5’-phosphate and 3’-hydroxyl group resulting from enzymatic cleavage by rna processing enzymes like dicer.
in the first step, rna adapters were ligated onto each end of the rna molecules, and a reverse transcription reaction was used to create single-stranded cdna.
this cdna was then pcr amplified with a universal primer and a second primer containing one of  <dig> uniquely indexed tags to allow multiplexing.
size selection of the cdna constructs was performed using a 3 % gel cassette on the pippin prep  to reduce the library to mature mirnas and other regulatory rnas in the 20–30 bp size range and to remove adapter-adapter products.
the resulting cdna libraries then underwent a quality check on the bioanalyzer hs dna assay  to confirm the final library size and on the agilent mx3005p quantitative pcr machine using the kapa library quantification kit  to determine concentration.
a  <dig> nm stock was created, and samples were pooled by molarity for multiplexing.
from the pool,  <dig> pm was loaded into each well for the flow cell on the illumina cbot for cluster generation.
the flow cell was then loaded onto the illumina hiseq  <dig> utilizing v <dig> chemistry and hta  <dig> .
the raw sequencing reads in bcl format were processed through casava- <dig> .
<dig> for fastq conversion and demultiplexing.
the rta chastity filter was used, and only the pf  reads were retained for further analysis.
sequencing data analysis
we implemented a custom analysis pipeline for small rna sequencing data processing.
we used cutadapt  <cit>  to trim 3’ adapters for raw reads.
quality control on raw data was performed using qc <dig>  <cit> .
all reads with length less than  <dig> were discarded.
the adaptor-trimmed reads were formatted into a non-redundant fastq file where the read sequence and copy number was recorded for each unique tag.
the usable unique reads were mapped to the whole genome using bowtie <dig>  <cit>  with only one mismatch allowed.
the latest sanger microrna database, mirbase <dig>  <cit> , was used to quantify microrna isomirs by reads mapped with position  <dig>  + <dig>  and + <dig> from the 3’ terminal of the mirbase locus.
the trna database was prepared by combining the latest ucsc trna database with the trna loci of mitochondria from the ensembl database  <cit> .
the reads mapped with trna loci were used not only for trna quantification, but also for trna mapping position coverage analysis.
the reads mapped to mirnas, trnas, and other small rnas  were used for response category analysis.
in addition, single nucleotide variations  can also be detected through rna-seq methods .
we identified snvs in the trna reads using varscan  <dig>  <cit>  and filtered the snvs based on the varscan  <dig> recommended filters.
inconsistent snvs between paired samples were removed from analysis.
only snvs with zero appearance in control samples were included in the analysis.
response analysis
a linear regression model was fit using response scores as outcomes and trna expression as predictors.
differential expression between mds biopsy samples and control samples were conducted using multirankseq  <cit> , which is based on the combination of results from three distinct rna-seq analysis packages: deseq  <cit> , edger  <cit> , and bayseq  <cit> .
a trna was considered significantly differentially expressed if the adjusted p-values for all three analysis packages were less than  <dig> .
the elastic net method  <cit>  was used to select a panel of trnas that together explain a large proportion of the variation in patient response.
elastic net regularization  <cit>  is a technique that conbines the l <dig> and l <dig> penalties of the lasso and ridge regression methods.
the elastic net method provides variable selection to produce parsimonious and interpretable models in the p > > n case without being severely limited by the sample size.
this method greatly improves performance in the case of highly correlated predictor variables , through the identification of groups of phenotypes with significantly high correlation that contribute the most to the variation in the data.
unsupervised hierarchical cluster analysis was conducted using heatmap <dig>  <cit> .
functional analysis of trnas/tdrs was carried out manually by searching online due to the lack of a functional database for trna.
in addition, a fisher’s exact test was used to identify snvs associated with mds.
RESULTS
we found that for the reads obtained using a standard mirna-seq protocol, the number of reads aligned to tdrs  vastly outnumbered those aligning to mirnas  .
the read length distribution, after trimming the adaptor with cutadapt  <cit> , shows a strong modal peak at  <dig> nucleotides, with the expected  <dig> nucleotide mean length of mirnas  forming only a small secondary peak in the distribution.
the peak at  <dig> bases indicates the abundance of trna species that have been cleaved in the anticodon loop from their full length of 73– <dig> nucleotides  <cit> .
the complete sequenced position of each trna is provided in additional file 2: figure s <dig>  numerous small rna mapping statistics can be located in  .fig.
<dig> read count and alignment distribution example taken from one sample.
the figures were produced using all read counts per category, not just unique reads per category.
the other samples in this study follow a similar pattern.
a. read count distribution after trimming adaptors.
the smaller peak at  <dig> base pairs indicates the abundance of mirna and the larger peak at  <dig> base pairs indicates the abundance primarily of trna.
b. the reads alignment distribution by rna type.
the majority of the reads aligned to trna instead of mirna
using multirankseq, we identified six tdrs that are significantly differentially expressed between mds samples and control samples .
three tdrs demonstrated increased expression in mds  while three were decreased tca tca), and chr <dig> trna4-thragt).
the complete results can be viewed in additional file 1: table s <dig>  unsupervised cluster analysis showed that trna profiling was able to distinguish between mds and control samples .
no tdrs were significantly differentially expressed between pre- and post-treatment samples .table  <dig> differentially expressed trna derivatives
alog2fc = log  <dig> fold change
bpadj = adjusted p -value
cdeseq <dig>  edger, and bayseq are the three rnaseq differential expression analysis packages used in this analysis.
bayseq does not generate fold change, thus no fold change from bayseq can be reported
fig.
<dig>
a. cluster analysis and heatmap using trna expression of all samples.
three phenotype bars are drawn below the dendrogram: pre-treatment, post-treatment and normal controls.
two clusters are visible .
these two clusters do not separate pre- and post-treatment, but distinguish mds and normal samples reasonably well.
b. the six differentially expressed trna between disease and normal
using the elastic net with ten-fold cross-validation  for selecting the model with the smallest cv error, we identified a panel of four trna fragments  whose combined expression in the pre-treatment samples together are predictive of the likelihood of response.
by fitting a multivariable linear regression model using expression values of these four tdrs, two out of the four trna species showed significant associations with response  adjusted for other variables in the model.
overall, the model explained roughly 67 % of the variation in treatment response  .
it should be noted that the first  <dig> nucleotides of chr <dig> trna <dig> valtac and chrx.trna <dig> valtac show perfect sequence homology, so these different trnas cannot be distinguished by these methods.
therefore, it is possible that one or the other species are predictive of response rather than both.table  <dig> linear regression results on trnas and treatment association
notice that not all trna are significant in this table; however, acting together, they explain the greatest amount of variation in treatment
aeffect is explained as per  <dig> unit increase in trna expression; the treatment changes the effect amount of unit
using the change of tdr expression between pre- and post-treatment , we conducted the same statistical analysis.
three tdrs  were identified as predictors for treatment response.
multivariate linear regression was fit using treatment response as the outcome and tdr expression changes in the three tdrs as predictors.
pseudottt trended towards an association with treatment response .
roughly 40 % of the variation in patient responses can be explained using the difference of expression between pre- and post-treatment of these three tdrs  .
it should be noted that the first  <dig> nucleotides of chr <dig> trna <dig> glygcc and chr <dig> trna <dig>  glygcc are also homologous, leading to alignment ambiguity.
functional correlation of trna expression is challenging due to the dearth of available public resources focused on trnas.
we focused on mitochondrial trnas  as a limited subset of trna species to study, in particular those mt-trnas that have a demonstrated disease association in the literature .
in table  <dig>  we show the false discovery rate  p-values for the ability of these tdrs from differential expression analysis methods to distinguish between mds and control  and their association with treatment response  for all  <dig> mt- tdrs.
thirteen of the  <dig> mt-tdrs had significantly adjusted p-values from at least two of the three methods from multirankseq in the discrimination of mds from controls.
mt-ts <dig> was the only trna from this subset that had significant association with treatment response.table  <dig> test statistics of previously identified mitochondria trnas with disease associations
through snv analysis, we identified  <dig> snvs that have potential associations with mds .
these  <dig> snvs were found in at least 20 % of mds samples and were not observed in any of the control samples.
only one out of the  <dig> snvs has been previously identified  by the  <dig> genomes project.
some of the snvs were identified in the same trna with the same anticodon but were located at different genomic positions.
it is possible that they are the same snv duplicated by alignment ambiguity.table  <dig> trna single nucleotide variant analysis
discussion
in this study, using mirna sequencing data generated from  <dig> bm-mnc samples and a novel bioinformatics approach, we demonstrated that tdrs are a major byproduct of mirna sequencing, exceeding the abundance of mirnas detected by this methodology by nearly 18-fold.
this finding demonstrates the difficulties inherent in studying low-abundance rna species by expression analysis and emphasizes the importance of adequate filtration.
due to the potential information embodied in the tdr dataset and a lack of published data on the topic, we chose to study the expression of tdrs in mds.
we demonstrated that six tdrs are differentially found in mds and control samples, and that a pattern of trna species could be used to differentiate mds and control samples using unsupervised hierarchical clustering.
most significantly, the expression of tdrs in pre-treatment samples was found to predict response to treatment with dnmtis.
these tdrs were found to explain 67 % of the variation in treatment response.
the expression of these entities could be used in the clinical setting to select patients likely to respond to dnmtis.
since only a minority of patients achieve a clinically meaningful response of hematologic improvement, partial response, or complete response to dnmti therapy, and since it typically takes 3– <dig> cycles before a response is evident  <cit> , pre-treatment identification of patients likely not to benefit from dnmti therapy would enable earlier decisions about stem cell transplantation or investigational therapies in those patients.
using changes in tdr expression between pre- and post-treatment samples in model fitting would explain less variation in the data ; however, since pre-treatment trna expression is a more clinically useful analysis, the better prediction value is fortuitous.
unfortunately, due to the small size of this study, internal validation methods have limited utility.
the optimal method to demonstrate the sensitivity and specificity of these biomarkers is to utilize an independent dataset.
this will be the focus of future efforts.
due to the recent developments in the study of mt-trnas and their association with various human diseases  <cit> , the intriguing studies on mitochondrial aberrations in mds  <cit> , and the identification of chrm.trna <dig> ts <dig> and chrm.trna <dig> tc in our initial analyses, we examined the expression of mt-trna-derived species in mds.
we identified  <dig> mt-tdrs which individually are differentially expressed in mds compared to control samples, as well as in responders and non-responders to dnmti therapy.
again, these may be valuable diagnostic and prognostic markers in mds; however, the correlation between this functional subset tdrs and the mitochondrial aberrations seen in mds is unclear.
finally, we found increased snvs within the trna sequences in the mds samples as compared to the controls.
while this may simply point to underlying genomic instability of myeloid neoplasms, studies in acute myeloid leukemia have not generally found increased absolute numbers of mutations  <cit> , and many mutations are found as a result of aging alone in individuals without any significant hematolymphoid malignancy .
therefore, there may be real functional consequences of these snvs which may affect the functions of these trna-derived species.
it should be noted that rna-seq is not a well-accepted method for the identification of snvs due to the high error rate of the reverse transcription step and the rna editing process  <cit> .
in addition, the presence of numerous modified nucleotides in trnas  can result in base errors during the reverse transcription step of rna-seq  <cit> .
however, these snvs were not identified to an appreciable level in the control samples, which were processed identically in a blinded fashion, suggesting that the significance of this finding is not an artifact.
study of trna fragments using mirna sequencing data is not without limitation.
as we described previously, trnas are typically 73– <dig> nucleotides in length, composed of three hairpin-turns in a cloverleaf two-dimensional structure.
the most 5’ of these loops is the d loop.
within this loop there is a site that has been shown to be sensitive to cleavage either by dicer or by simple incubation with mgcl <dig> to afford a  <dig> nucleotide product from the 5’ end  <cit> .
however, figures in this same publication also identify a cleavage product of greater than  <dig> nucleotides  <cit>  that is too small to correspond to the residual 3’ fragment.
as we have shown,  <dig> nucleotide trna derivatives were identified by our sequencing efforts, resulting from cleavage within the second hairpin loop containing the anticodon.
our size selection step would have excluded most of the  <dig> nucleotide product.
in addition, cleavage in the anti-codon loop occurs physiologically as a result of cleavage by angiogenin  <cit>  and may be targeted by other nucleases as well  <cit> .
lastly, due to the poor reverse transcription of full length parent trnas due to nucleotide modifications, apparent tdrs may result from incomplete reverse transcription  <cit> .
thus, for various biological and technical reasons, it is possible that under the conditions of this study, the  <dig> base derivative may be the predominant 5’ tdr product.
the lack of knowledge of the entire sequence of the trna creates ambiguous trna annotation.
many of the trna isotypes have not only the same or similar anticodons, but also highly homologous sequences 5’ to the anticodon.
this creates challenges in the identification of specific species that may be diagnostically or prognostically useful.
as mentioned in the results, several of the prognostic trna fragments are homologous.
therefore, although the degenerate species are individually identified, the prognostic power may lie in only one of the entities.
also, in table  <dig>  a g > a snv was identified in five seraga trnas.
since the first  <dig> base pairs of these five trnas are the same, the redundant alignment resulted in the identification of a g > a snv in all five seraga trnas for the homologous position, when in fact it may have only been seen in as few as one of the five seraga isotypes.
the actual nucleotide sequences of each of the tdrs specifically mentioned in this study are available for comparison in additional file 1: table s <dig>  furthermore, we used varscan  <dig> to identify snvs.
typically, at least 20 % of the reads need to support the alternative allele in order to be called an snv by varscan  <dig>  so, although none of the control samples in table  <dig> contained an snv as determined by varscan  <dig>  there may have been low levels of reads that did not meet the 20 % threshold for identification of an snv.
it is unclear if the expression of tdrs represents the true expression of the full-length trnas themselves, if they are physiologic by-products whose presence at different levels suggests differential processing and, therefore, half-lives of the trnas, or if they represent biologically active entities in themselves.
needless to say, this ambiguity complicates the further exploration of these results, requiring novel methods to unravel this unknown biology.
recently, new approaches to library preparation have emerged, including the use of demethylases and other strategies to remove modifications prior to reverse transcription.
moreover, the use of group ii intron reverse transcriptases, that bind more tightly than retroviral enzymes to the rna template, has been demonstrated to overcome modification barriers to afford full-length trna cdnas in rna-seq libraries .
these methods may be useful in the future to distinguish the biological tdrs from those that result from technical artifact.
the exact mechanism by which tdrs may play a role in the diagnosis of mds and prognostication of response to dnmtis is unclear.
there are multitudinous effects of global hypomethylation that results from dnmti therapy.
some tdrs have been demonstrated to function independently as biologically active entities.
lee et al.
have identified small trna-derived species, including trf- <dig>  which promotes the g2/m transition  <cit> .
however, the species identified by lee et al.
are typically 17– <dig> bases long, smaller than the  <dig> nucleotide reads identified in this study.
several reports suggest that 30– <dig> nucleotide long trna-derived fragments may play a role in the biogenesis of other small rnas, and these have been found in wide ranges of species, including humans.
yamasaki et al.
suggest that these fragments may mediate stress-related translational repression  <cit> .
however, these explanations would have unclear biologic importance in the differential expression of specific trna-derived sequences.
alterations in the activity of dicer, and hence of mirnas, are a common feature in cancers .
trnas have been proposed as alternative dicer <dig> substrates  <cit> .
whether certain trnas are selectively targeted by dicer <dig> in cancer in general, and mds in particular, remains to be seen, but it provides an intriguing hypothesis for the differentially expressed tdrs identified in this study in association with both the diagnosis and prognosis of mds.
CONCLUSIONS
this study suggests that trna-derived fragment sequencing can provide an additional source of data that potentially provides clinically useful applications in diagnosis and prognosis of disease.
in addition, this study raises intriguing questions about the biology of these trna derivatives in mds and other cancers.
additional files
additional file 1: table s <dig>  sample description.
table s <dig>  phenotype and alignment results for all samples.
table s <dig>  differential expression analysis between mds and control samples.
table s <dig>  differential expression analysis between pre and post treatment.
table s <dig>  nucleotide sequences of the trnas described in this paper.
additional file 2: figure s <dig>  using a single sample as an example, this figure shows the alignment position of each trna.
the majority of the alignment lies between residues 1– <dig> from the 5’ end.
the size of each dot indicates the abundance of reads aligned to this position.
abbreviations
mirnamicrorna
mirna-seqmicrorna sequencing
mdsmyelodysplastic syndromes
trnatransfer rna
trhstransfer rna halves
tdrstransfer rna-derived small rnas
ntnucleotide
amlacute myeloid leukemia
dnmtisdna methyltransferase inhibitors
bm-mncsbone marrow mononuclear cells
totrnatotal rna
snvsingle nucleotide variant
cvcross-validation
mt-trnasmitochondrial transfer rnas
fdrfalse discovery rate
competing interests
the authors declare that they have no competing interests related to this research.
authors’ contributions
yg and ask contributed to the conception and design of the study.
ab, sm, be, mj, sas, eag, and ask were involved in the acquisition of the data.
yg, ab, sm, qs, fy, c-il, pfs, sz, and ask were involved in the analysis and interpretation of the data.
yg, ab, fy, kcv, and ask were involved in the drafting of the manuscript and critical revisions for important intellectual content.
all authors have given final approval of the version to be published.
yg and ask agree to be accountable for all aspects of the work.
authors’ information
not applicable.
availability of data and materials
not applicable.

 rna next generation sequencing technology has replaced microarrays as the primary platform for gene expression profiling. this same technology has also been used to study the expression of micrornas, although there are fewer direct comparisons between mirna arrays and mirna-seq. during mirna-seq library preparation, small rnas are selected by electrophoresis on a size typically ranging from to nucleotides. this range of size selection allows for the capture of many other species of small rnas in addition to mirnas. one of the major byproducts of size selection for mirnas is trna-derived small rnas. this may be the result of either active cleavage or an artifact of the mirna-seq library construction. full length or parent trnas are adaptor molecules composed of rna with a length typically ranging from to nts. they serve as the physical link that translates messenger rna into protein. cleavage of trnas by the rnase iii enzyme, angiogenin, may occur in a number of reactive conditions to produce trna-derived halves. likewise, trnas can also be cleaved in a dicer-dependent manner or as an in vitro phenomenon by incubation with mgcl or nuclease s. all of these processes result in cleavage within a hairpin loop. d loop cleavage results in a 19 nt trna-derived fragment, and cleavage of the anticodon loop before the anti-codon affords a 33 nt trh. these 5’ tdr products, in particular the 33-base trh, may be captured by mirna-seq. in addition, apparent tdrs and trnas with base errors found in rna-seq datasets may arise from the extensive chemical modifications of parent trna nucleotides. while these modifications contribute physiologically to rna stability, function, and structure, they may be either misincorporated during the requisite reverse transcription step, resulting in specific base errors, or result in blockage of reverse transcription, with production of a truncated cdna product. therefore, sequencing datasets may contain many reads that appear to be tdrs, but in fact are actually derived from full-length parent trna templates. the quantification of trnas and tdrs warrant study through high-throughput sequencing. aberrations in trnas have been linked to a range of diseases, including neurological disorders, cancer, type ii diabetes, and mitochondrial diseases. in addition, tdrs have been identified in a number of cancer cell lines, and some tdrs may play defined biological roles, such as the role of trf- in proliferation. thus far, only one study has found a link between somatic mutations in mitochondrial trna, aka mt-tl1) and myelodysplastic syndromes. the potential functionality of trna in mds remains largely unexplored. as part of our efforts to study mds, we have previously examined mirna expression in mds. mds is a common disease of the elderly characterized by ineffective maturation of hematopoietic cells. this ineffective maturation manifests clinically as low blood counts and morphologically as dysplasia. the mds state has been associated with expression changes in mirnas. several entities have been discovered on 5q, which contribute to the disease phenotype in the 5q minus syndrome subtype of mds. mir- and mir- have recently been identified in the commonly deleted region, and their knockdown appears to be associated with the thrombocytosis and dysmegakaryopoiesis seen in 5q minus syndrome. researchers, including those in our group, have identified several other mirnas whose expression is dysregulated with the diagnosis of other subtypes of mds. published research on other types of small rnas in mds is limited. mds has a median survival of only 18– months, with death resulting from either cytopenic complications or transformation to acute myeloid leukemia. currently, there are only two fda-approved drugs for the treatment of all types of mds, both of which are dna methyltransferase inhibitors: 5-azacytidine, and decitabine. however, only 40–47 % of patients achieve the clinically meaningful responses of hematologic improvement or better with dnmtis. although several recent studies have suggested the use of mutational analysis by dna next generation sequencing, genetics, or the presence or absence of key pharmacodynamic and pharmacokinetic markers to predict response, as yet, there are no clinical or laboratory parameters in routine clinical practice that accurately predict response to dnmtis. the methylation status of specific loci has been shown to be prognostically relevant in the treatment of a closely related disease, chronic myelomonocytic leukemia, by decitabine. in this study, we examined the expression of trnas in paired pre- and post-treatment samples from patients with the diagnosis of mds on receiving treatment with dnmtis. in mirna sequencing data, the abundance of trna fragments vastly outnumbers mirnas. this gives us a great opportunity to study trna expression using mirna sequencing data through their anticodon sequences. methods ethics statement this study was conducted in accordance with the declaration of helsinki with the approval of the institutional review boards at roswell park cancer institute and vanderbilt university medical center. consent statement written informed consent for the patient-derived specimens was obtained prior to the acquisition in all cases. the consents were approved by the institutional review boards at roswell park cancer institute and vanderbilt university medical center. specimen description bone marrow aspirate specimens were obtained from the roswell park cancer institute with appropriate approval by the institutional review boards at roswell park and vanderbilt university medical center. twenty-two of the bone marrow samples were from control patients. the other samples were from mds patients. out of mds samples, were pre-treatment specimens, and were obtained after treatment with 5-azacytidine or decitabine, both dnmtis. for each of the pre- and post-treatment paired samples, a response score was obtained based on clinical and pathological criteria to indicate the effectiveness of treatment. bone marrow mononuclear cells were isolated from fresh, unsorted bone marrow aspirate specimens by the ficoll method using cellgro lymphocyte separation medium. after performing a cell count, the cells were re-pelleted from dulbecco’s phosphate buffered saline solution and resuspended in gibco recovery freezing medium, and then frozen gradually, prior to storage in liquid nitrogen. cells were frozen at a density of to million cells per ml. rna isolation and rna sequencing total rna was isolated using a mirvana mirna isolation kit per manufacturer’s instructions. sequencing library construction was performed on the totrna from all samples, each obtained from a single bone marrow aspirate specimen. rnas were isolated by the mirvana rna isolation kit, which was found to perform better than the trizol mirna isolation kit. libraries were prepared using the truseq small rna sample preparation kit. the small rna protocol specifically ligates rna adapters to mature mirnas that have a 5’-phosphate and 3’-hydroxyl group resulting from enzymatic cleavage by rna processing enzymes like dicer. in the first step, rna adapters were ligated onto each end of the rna molecules, and a reverse transcription reaction was used to create single-stranded cdna. this cdna was then pcr amplified with a universal primer and a second primer containing one of uniquely indexed tags to allow multiplexing. size selection of the cdna constructs was performed using a 3 % gel cassette on the pippin prep to reduce the library to mature mirnas and other regulatory rnas in the 20–30 bp size range and to remove adapter-adapter products. the resulting cdna libraries then underwent a quality check on the bioanalyzer hs dna assay to confirm the final library size and on the agilent mx3005p quantitative pcr machine using the kapa library quantification kit to determine concentration. a nm stock was created, and samples were pooled by molarity for multiplexing. from the pool, pm was loaded into each well for the flow cell on the illumina cbot for cluster generation. the flow cell was then loaded onto the illumina hiseq utilizing v chemistry and hta. the raw sequencing reads in bcl format were processed through casava-. for fastq conversion and demultiplexing. the rta chastity filter was used, and only the pf reads were retained for further analysis. sequencing data analysis we implemented a custom analysis pipeline for small rna sequencing data processing. we used cutadapt to trim 3’ adapters for raw reads. quality control on raw data was performed using qc. all reads with length less than were discarded. the adaptor-trimmed reads were formatted into a non-redundant fastq file where the read sequence and copy number was recorded for each unique tag. the usable unique reads were mapped to the whole genome using bowtie with only one mismatch allowed. the latest sanger microrna database, mirbase, was used to quantify microrna isomirs by reads mapped with position + and + from the 3’ terminal of the mirbase locus. the trna database was prepared by combining the latest ucsc trna database with the trna loci of mitochondria from the ensembl database. the reads mapped with trna loci were used not only for trna quantification, but also for trna mapping position coverage analysis. the reads mapped to mirnas, trnas, and other small rnas were used for response category analysis. in addition, single nucleotide variations can also be detected through rna-seq methods. we identified snvs in the trna reads using varscan and filtered the snvs based on the varscan recommended filters. inconsistent snvs between paired samples were removed from analysis. only snvs with zero appearance in control samples were included in the analysis. response analysis a linear regression model was fit using response scores as outcomes and trna expression as predictors. differential expression between mds biopsy samples and control samples were conducted using multirankseq, which is based on the combination of results from three distinct rna-seq analysis packages: deseq, edger, and bayseq. a trna was considered significantly differentially expressed if the adjusted p-values for all three analysis packages were less than. the elastic net method was used to select a panel of trnas that together explain a large proportion of the variation in patient response. elastic net regularization is a technique that conbines the l and l penalties of the lasso and ridge regression methods. the elastic net method provides variable selection to produce parsimonious and interpretable models in the p > > n case without being severely limited by the sample size. this method greatly improves performance in the case of highly correlated predictor variables, through the identification of groups of phenotypes with significantly high correlation that contribute the most to the variation in the data. unsupervised hierarchical cluster analysis was conducted using heatmap. functional analysis of trnas/tdrs was carried out manually by searching online due to the lack of a functional database for trna. in addition, a fisher’s exact test was used to identify snvs associated with mds.  we found that for the reads obtained using a standard mirna-seq protocol, the number of reads aligned to tdrs vastly outnumbered those aligning to mirnas. the read length distribution, after trimming the adaptor with cutadapt, shows a strong modal peak at nucleotides, with the expected nucleotide mean length of mirnas forming only a small secondary peak in the distribution. the peak at bases indicates the abundance of trna species that have been cleaved in the anticodon loop from their full length of 73– nucleotides. the complete sequenced position of each trna is provided in additional file 2: figure s numerous small rna mapping statistics can be located in.fig. read count and alignment distribution example taken from one sample. the figures were produced using all read counts per category, not just unique reads per category. the other samples in this study follow a similar pattern. a. read count distribution after trimming adaptors. the smaller peak at base pairs indicates the abundance of mirna and the larger peak at base pairs indicates the abundance primarily of trna. b. the reads alignment distribution by rna type. the majority of the reads aligned to trna instead of mirna using multirankseq, we identified six tdrs that are significantly differentially expressed between mds samples and control samples. three tdrs demonstrated increased expression in mds while three were decreased tca tca), and chr trna4-thragt). the complete results can be viewed in additional file 1: table s unsupervised cluster analysis showed that trna profiling was able to distinguish between mds and control samples. no tdrs were significantly differentially expressed between pre- and post-treatment samples.table differentially expressed trna derivatives alog2fc = log fold change bpadj = adjusted p -value cdeseq edger, and bayseq are the three rnaseq differential expression analysis packages used in this analysis. bayseq does not generate fold change, thus no fold change from bayseq can be reported fig. a. cluster analysis and heatmap using trna expression of all samples. three phenotype bars are drawn below the dendrogram: pre-treatment, post-treatment and normal controls. two clusters are visible. these two clusters do not separate pre- and post-treatment, but distinguish mds and normal samples reasonably well. b. the six differentially expressed trna between disease and normal using the elastic net with ten-fold cross-validation for selecting the model with the smallest cv error, we identified a panel of four trna fragments whose combined expression in the pre-treatment samples together are predictive of the likelihood of response. by fitting a multivariable linear regression model using expression values of these four tdrs, two out of the four trna species showed significant associations with response adjusted for other variables in the model. overall, the model explained roughly 67 % of the variation in treatment response. it should be noted that the first nucleotides of chr trna valtac and chrx.trna valtac show perfect sequence homology, so these different trnas cannot be distinguished by these methods. therefore, it is possible that one or the other species are predictive of response rather than both.table linear regression results on trnas and treatment association notice that not all trna are significant in this table; however, acting together, they explain the greatest amount of variation in treatment aeffect is explained as per unit increase in trna expression; the treatment changes the effect amount of unit using the change of tdr expression between pre- and post-treatment, we conducted the same statistical analysis. three tdrs were identified as predictors for treatment response. multivariate linear regression was fit using treatment response as the outcome and tdr expression changes in the three tdrs as predictors. pseudottt trended towards an association with treatment response. roughly 40 % of the variation in patient responses can be explained using the difference of expression between pre- and post-treatment of these three tdrs. it should be noted that the first nucleotides of chr trna glygcc and chr trna glygcc are also homologous, leading to alignment ambiguity. functional correlation of trna expression is challenging due to the dearth of available public resources focused on trnas. we focused on mitochondrial trnas as a limited subset of trna species to study, in particular those mt-trnas that have a demonstrated disease association in the literature. in table  we show the false discovery rate p-values for the ability of these tdrs from differential expression analysis methods to distinguish between mds and control and their association with treatment response for all mt- tdrs. thirteen of the mt-tdrs had significantly adjusted p-values from at least two of the three methods from multirankseq in the discrimination of mds from controls. mt-ts was the only trna from this subset that had significant association with treatment response.table test statistics of previously identified mitochondria trnas with disease associations through snv analysis, we identified snvs that have potential associations with mds. these snvs were found in at least 20 % of mds samples and were not observed in any of the control samples. only one out of the snvs has been previously identified by the genomes project. some of the snvs were identified in the same trna with the same anticodon but were located at different genomic positions. it is possible that they are the same snv duplicated by alignment ambiguity.table trna single nucleotide variant analysis discussion in this study, using mirna sequencing data generated from bm-mnc samples and a novel bioinformatics approach, we demonstrated that tdrs are a major byproduct of mirna sequencing, exceeding the abundance of mirnas detected by this methodology by nearly 18-fold. this finding demonstrates the difficulties inherent in studying low-abundance rna species by expression analysis and emphasizes the importance of adequate filtration. due to the potential information embodied in the tdr dataset and a lack of published data on the topic, we chose to study the expression of tdrs in mds. we demonstrated that six tdrs are differentially found in mds and control samples, and that a pattern of trna species could be used to differentiate mds and control samples using unsupervised hierarchical clustering. most significantly, the expression of tdrs in pre-treatment samples was found to predict response to treatment with dnmtis. these tdrs were found to explain 67 % of the variation in treatment response. the expression of these entities could be used in the clinical setting to select patients likely to respond to dnmtis. since only a minority of patients achieve a clinically meaningful response of hematologic improvement, partial response, or complete response to dnmti therapy, and since it typically takes 3– cycles before a response is evident, pre-treatment identification of patients likely not to benefit from dnmti therapy would enable earlier decisions about stem cell transplantation or investigational therapies in those patients. using changes in tdr expression between pre- and post-treatment samples in model fitting would explain less variation in the data; however, since pre-treatment trna expression is a more clinically useful analysis, the better prediction value is fortuitous. unfortunately, due to the small size of this study, internal validation methods have limited utility. the optimal method to demonstrate the sensitivity and specificity of these biomarkers is to utilize an independent dataset. this will be the focus of future efforts. due to the recent developments in the study of mt-trnas and their association with various human diseases, the intriguing studies on mitochondrial aberrations in mds, and the identification of chrm.trna ts and chrm.trna tc in our initial analyses, we examined the expression of mt-trna-derived species in mds. we identified mt-tdrs which individually are differentially expressed in mds compared to control samples, as well as in responders and non-responders to dnmti therapy. again, these may be valuable diagnostic and prognostic markers in mds; however, the correlation between this functional subset tdrs and the mitochondrial aberrations seen in mds is unclear. finally, we found increased snvs within the trna sequences in the mds samples as compared to the controls. while this may simply point to underlying genomic instability of myeloid neoplasms, studies in acute myeloid leukemia have not generally found increased absolute numbers of mutations, and many mutations are found as a result of aging alone in individuals without any significant hematolymphoid malignancy. therefore, there may be real functional consequences of these snvs which may affect the functions of these trna-derived species. it should be noted that rna-seq is not a well-accepted method for the identification of snvs due to the high error rate of the reverse transcription step and the rna editing process. in addition, the presence of numerous modified nucleotides in trnas can result in base errors during the reverse transcription step of rna-seq. however, these snvs were not identified to an appreciable level in the control samples, which were processed identically in a blinded fashion, suggesting that the significance of this finding is not an artifact. study of trna fragments using mirna sequencing data is not without limitation. as we described previously, trnas are typically 73– nucleotides in length, composed of three hairpin-turns in a cloverleaf two-dimensional structure. the most 5’ of these loops is the d loop. within this loop there is a site that has been shown to be sensitive to cleavage either by dicer or by simple incubation with mgcl to afford a nucleotide product from the 5’ end. however, figures in this same publication also identify a cleavage product of greater than nucleotides that is too small to correspond to the residual 3’ fragment. as we have shown, nucleotide trna derivatives were identified by our sequencing efforts, resulting from cleavage within the second hairpin loop containing the anticodon. our size selection step would have excluded most of the nucleotide product. in addition, cleavage in the anti-codon loop occurs physiologically as a result of cleavage by angiogenin and may be targeted by other nucleases as well. lastly, due to the poor reverse transcription of full length parent trnas due to nucleotide modifications, apparent tdrs may result from incomplete reverse transcription. thus, for various biological and technical reasons, it is possible that under the conditions of this study, the base derivative may be the predominant 5’ tdr product. the lack of knowledge of the entire sequence of the trna creates ambiguous trna annotation. many of the trna isotypes have not only the same or similar anticodons, but also highly homologous sequences 5’ to the anticodon. this creates challenges in the identification of specific species that may be diagnostically or prognostically useful. as mentioned in the results, several of the prognostic trna fragments are homologous. therefore, although the degenerate species are individually identified, the prognostic power may lie in only one of the entities. also, in table  a g > a snv was identified in five seraga trnas. since the first base pairs of these five trnas are the same, the redundant alignment resulted in the identification of a g > a snv in all five seraga trnas for the homologous position, when in fact it may have only been seen in as few as one of the five seraga isotypes. the actual nucleotide sequences of each of the tdrs specifically mentioned in this study are available for comparison in additional file 1: table s furthermore, we used varscan to identify snvs. typically, at least 20 % of the reads need to support the alternative allele in order to be called an snv by varscan so, although none of the control samples in table  contained an snv as determined by varscan there may have been low levels of reads that did not meet the 20 % threshold for identification of an snv. it is unclear if the expression of tdrs represents the true expression of the full-length trnas themselves, if they are physiologic by-products whose presence at different levels suggests differential processing and, therefore, half-lives of the trnas, or if they represent biologically active entities in themselves. needless to say, this ambiguity complicates the further exploration of these results, requiring novel methods to unravel this unknown biology. recently, new approaches to library preparation have emerged, including the use of demethylases and other strategies to remove modifications prior to reverse transcription. moreover, the use of group ii intron reverse transcriptases, that bind more tightly than retroviral enzymes to the rna template, has been demonstrated to overcome modification barriers to afford full-length trna cdnas in rna-seq libraries. these methods may be useful in the future to distinguish the biological tdrs from those that result from technical artifact. the exact mechanism by which tdrs may play a role in the diagnosis of mds and prognostication of response to dnmtis is unclear. there are multitudinous effects of global hypomethylation that results from dnmti therapy. some tdrs have been demonstrated to function independently as biologically active entities. lee et al. have identified small trna-derived species, including trf- which promotes the g2/m transition. however, the species identified by lee et al. are typically 17– bases long, smaller than the nucleotide reads identified in this study. several reports suggest that 30– nucleotide long trna-derived fragments may play a role in the biogenesis of other small rnas, and these have been found in wide ranges of species, including humans. yamasaki et al. suggest that these fragments may mediate stress-related translational repression. however, these explanations would have unclear biologic importance in the differential expression of specific trna-derived sequences. alterations in the activity of dicer, and hence of mirnas, are a common feature in cancers. trnas have been proposed as alternative dicer substrates. whether certain trnas are selectively targeted by dicer in cancer in general, and mds in particular, remains to be seen, but it provides an intriguing hypothesis for the differentially expressed tdrs identified in this study in association with both the diagnosis and prognosis of mds.  this study suggests that trna-derived fragment sequencing can provide an additional source of data that potentially provides clinically useful applications in diagnosis and prognosis of disease. in addition, this study raises intriguing questions about the biology of these trna derivatives in mds and other cancers. additional files additional file 1: table s sample description. table s phenotype and alignment results for all samples. table s differential expression analysis between mds and control samples. table s differential expression analysis between pre and post treatment. table s nucleotide sequences of the trnas described in this paper. additional file 2: figure s using a single sample as an example, this figure shows the alignment position of each trna. the majority of the alignment lies between residues 1– from the 5’ end. the size of each dot indicates the abundance of reads aligned to this position. abbreviations mirnamicrorna mirna-seqmicrorna sequencing mdsmyelodysplastic syndromes trnatransfer rna trhstransfer rna halves tdrstransfer rna-derived small rnas ntnucleotide amlacute myeloid leukemia dnmtisdna methyltransferase inhibitors bm-mncsbone marrow mononuclear cells totrnatotal rna snvsingle nucleotide variant cvcross-validation mt-trnasmitochondrial transfer rnas fdrfalse discovery rate competing interests the authors declare that they have no competing interests related to this research. authors’ contributions yg and ask contributed to the conception and design of the study. ab, sm, be, mj, sas, eag, and ask were involved in the acquisition of the data. yg, ab, sm, qs, fy, c-il, pfs, sz, and ask were involved in the analysis and interpretation of the data. yg, ab, fy, kcv, and ask were involved in the drafting of the manuscript and critical revisions for important intellectual content. all authors have given final approval of the version to be published. yg and ask agree to be accountable for all aspects of the work. authors’ information not applicable. availability of data and materials not applicable. 
the species of small rnas were enumerated, single nucleotide variants identified, and finally the differential expression of trna-derived species in the specimens correlated with diseasestatus and response to therapy. total rna was extracted from these specimens and submitted for next generation sequencing after an additional size exclusion step to enrich for small rnas. consequently, the libraries obtained for mirna sequencing also contain a variety of additional species of small rnas. this study looks at the prevalence of these other species obtained from bone marrow aspirate specimens and explores the predictive value of these small rnas in the determination of response to therapy in myelodysplastic syndromes. using mirna sequencing data generated from bone marrow aspirate samples of patients with known mds and controls, we demonstrated that transfer rna fragments are one of the most common species of small rna isolated from size selection. using the elastic net method, we identified four trnas-derived small rnas that together can explain 67 % of the variation in treatment response for mds patients. this study highlights the prevalence of tdrs in rna-seq studies focused on small rnas. paired pre and post treatment bone marrow aspirate specimens were obtained from patients with mds who were treated with either azacytidine or decitabine with additional non-mds control specimens. our analysis results suggest that trna fragments can be accurately detected through mirna sequencing data and that the expression of these species may be useful in the diagnosis of mds and the prediction of response to therapy. using trna expression values extracted from mirna sequencing data, we identified six trna fragments that are differentially expressed between mds and normal samples. however, as we have previously described, mirnas represent only a subset of the species obtained by size selection. finally, snvs within the tdrs were found in at least % of the mds samples and were not observed in any of the control specimens. although advances in sequencing technologies have popularized the use of microrna sequencing for the quantification of mirna expression, questions remain concerning the optimal methodologies for analysis and utilization of the data.
 although advances in sequencing technologies have popularized the use of microrna sequencing for the quantification of mirna expression, questions remain concerning the optimal methodologies for analysis and utilization of the data. the construction of a mirna sequencing library selects rna by length rather than type. however, as we have previously described, mirnas represent only a subset of the species obtained by size selection. consequently, the libraries obtained for mirna sequencing also contain a variety of additional species of small rnas. this study looks at the prevalence of these other species obtained from bone marrow aspirate specimens and explores the predictive value of these small rnas in the determination of response to therapy in myelodysplastic syndromes. methods paired pre and post treatment bone marrow aspirate specimens were obtained from patients with mds who were treated with either azacytidine or decitabine with additional non-mds control specimens. total rna was extracted from these specimens and submitted for next generation sequencing after an additional size exclusion step to enrich for small rnas. the species of small rnas were enumerated, single nucleotide variants identified, and finally the differential expression of trna-derived species in the specimens correlated with diseasestatus and response to therapy.  using mirna sequencing data generated from bone marrow aspirate samples of patients with known mds and controls, we demonstrated that transfer rna fragments are one of the most common species of small rna isolated from size selection. using trna expression values extracted from mirna sequencing data, we identified six trna fragments that are differentially expressed between mds and normal samples. using the elastic net method, we identified four trnas-derived small rnas that together can explain 67 % of the variation in treatment response for mds patients. similar analysis of specifically mitochondrial tdrs identified mt-tdrs which distinguished disease status in the samples and a single mt-tdr which predited response. finally, snvs within the tdrs were found in at least % of the mds samples and were not observed in any of the control specimens. discussion this study highlights the prevalence of tdrs in rna-seq studies focused on small rnas. the potential etiologies of these species, both technical and biologic, are discussed as well as important challenges in the interpretation of tdr data.  our analysis results suggest that trna fragments can be accurately detected through mirna sequencing data and that the expression of these species may be useful in the diagnosis of mds and the prediction of response to therapy. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keyword transfer rna microrna next generation sequencing predictive biomarkersdna methyltransferase inhibitorresponse to therapymyelodysplastic syndromes doris duke charitable foundation, clinical scientist development award112076kim annette s. national institutes of healthp ca068485guo yan issue-copyright-statement© the author 2015 