Analysis of microRNA transcriptome by deep sequencing of small RNA libraries of peripheral blood

1Small non-coding RNAs participate in a variety of processes from cell development and differentiation, stress responses to carcinogenesis by regulating gene expression. Regulatory non-coding RNAs have been reported from almost all organisms from bacteria to mammals. Among the various classes of non-coding small RNAs (sRNAs), the most conserved and prominent ones are the microRNAs or miRNAs [for a recent review see,]. Mature miRNA sequences are single stranded, typically 18-24 nucleotides long and encoded as a precursor molecule of about 60-120 nucleotides (in humans). These precursors are derived from processing pri-miRNA (usually in kilobases) by a ribonuclease, such as DROSHA. Pre-miRNAs are also further cleaved to generate active mature miRNAs with the help of DICER. So far more than 800 miRNAs have been described in human. miRNAs interact with 3'UTR of mRNAs through base pairing and bring about their degradation, destabilization, or repression of translation through RISC, a complex of multiple proteins and miRNA-mRNA adduct. Occasionally it can also up regulate gene expression.

2The expression profiles of miRNAs have been determined in order to understand the role of miRNAs in a specific biological process. The profiles generated revealed altered levels of miRNAs in different systems, such as oncogenesis and development. These studies showed strong correlation between specific miRNA expression and phenotype suggesting it can be used as a potential biomarker for diagnosis and prognosis in human cancer. Several groups of miRNAs have been identified that regulate the expression of tumor-associated genes while others seem to hold prognostic value in predicting patient survival. For example, miR-21 is frequently over expressed in various cancers. It is now considered as an oncomiR that acts by down regulating PTEN, a tumor suppressor gene. Similarly let-7 family of miRNAs is frequently down regulated in lung cancers. miR-15a and 16-1 are often deleted in chronic lymphocytic leukemia. These miRNAs target RAS, HMGA2 and BCL2 oncogenes, respectively, thereby regulating tumorigenesis. Expression of let-7a, miR-210 and miR-200 clusters was demonstrated to be a strong prognostic marker for lung cancer, breast cancer and ovarian cancer, respectively. Many of the fundamental processes are regulated by miRNAs, such as cell differentiation (miR-223, miR-145), apoptosis (miR-34, miR-16), body patterning (miR-9, miR-196), nervous system and muscle development (miR-134, miR-1 and miR-133).

3A number of different approaches have been used for expression profiling, such as northern analysis, cloning, real time polymerase chain reaction, microarray analysis and RNase protection assay. These methods are not generally useful for discovery and expression profiling of low-abundance transcripts or yet unidentified novel miRNAs. Recently a number of different platforms have been developed for carrying out large scale parallel sequencing in order to generate genome wide sequences in a short time and reduced cost. These have also been applied to analyze transcriptome including sRNA sequences. Generally sequences derived using these platforms can be processed to generate expression profiles of known genes and suitable computational methods can be employed to identify unknown genes. Some of these platforms along with custom computation pipelines have been used to study sRNAs from different systems, such as plants, human embryonic cells and developing chicken embryo. All these studies not only generated expression profiles of known miRNAs but also identified a few novel miRNAs. Morin et al. used this strategy in combination with RNAfold, MiPred and an in house SVM model to identify novel miRNAs. Though the study identified 104 novel miRNAs there was no attempt to verify these by other experiments. A software pipeline miRDeep has been developed to analyze large-scale sRNA sequencing data and identification of novel miRNAs. Some of the filters used in this approach are highly stringent and it is likely that many novel miRNAs may be missed. Generation of expression profiles and identification of novel miRNAs from deep sequencing is dependent on tools used for analysis. It appears from this discussion that there is a need for development of efficient automated pipelines for analysis of deep sequencing data and setting up validation pipeline for checking the results, as newer algorithms may allow improved profiling and identification of novel miRNAs.

4In this manuscript we have profiled sRNA expression from normal human peripheral blood mononuclear cells and two cancer cell lines using Solexa technology and developed automated computation pipelines for analyzing quantitative expression. Our pipelines use statistical analysis of the Solexa sequences for generation of expression profiles and a number of different methods for prediction of novel miRNAs. We believe that these pipelines are highly robust and can be useful for other studies.

5sRNAs were size selected by gel electrophoresis and then sequenced using a Solexa platform (Illumina, USA). The details of sequencing reads of the four samples are given in Additional file 1. The length of an average read was about 33-35 nucleotides. The processing of the raw sequences through computation pipeline is outlined in Figure 1 and described in "Methods". Briefly, the adaptor sequences were first removed from the sequence reads. Only those reads that were greater than 10 nucleotides were considered for further analysis. These sequences were clustered on the basis of sequence similarity and subjected to similarity searches using specific databases (rRNAs, tRNA, sn/snoRNAs, miRNAs, other non-coding RNAs). The reads that did not match known RNA sequences were checked to see if these were encoded by intergenic, intronic or exonic regions of the human genome (Figure 2).

6Sequencing reads derived from mRNAs and rRNAs made up about 2-4 and 2-13% of the total reads of all the samples, respectively except that from K562 cells. Since mRNAs and rRNAs are likely to have been derived by degradation it appears that the degradation was minimal in these samples. The sequencing reads derived from miRNAs were estimated to be about 60-80%, similar to some of the other studies on sRNAs. The result from K562 was quite different. Only about 18% of total sRNA population was found to be derived from miRNA. Low level of miRNAs in K562 cells is unlikely to be due to artifact as we got similar values in two independent experiments (data not shown). Out of a total of 904 miRNAs (718 major and 186 star sequences) present in miRBase version 14, 534 miRNAs were detected in at least one of the four samples sequenced. A total of 370 miRNAs were not detected in any of the samples.

7Absolute sequence reads were transformed into transcript abundance by first normalizing the data in 'transcripts per million (TPM)' for each library (see "Methods"). The expression levels ranged from less than 10 to more than 100,000 counts (Figure 3, refer Additional file 1). Thus the sequencing data revealed a wide range of expression levels spanning five orders of magnitude. Several members of the let-7 family, 103, 185 and 320a were observed as some of the highly expressed miRNAs (Figure 4). These miRNAs were also reported as highly abundant miRNAs in peripheral blood mononuclear cells using "Taqman microRNA assay". While the members of let-7 family, let-7a, f and g represented about 77% of the total miRNA counts, a number of miRNAs with less than 10 counts including singletons were also noticed (Additional file 1). It is not clear if the singletons represent true transcript or noise of the system. miR-219-5p was one of the singletons found in all the four samples suggesting that such sequences may not be due to experimental noise or a chance event. Star sequences of many miRNAs (20% of all known miRNAs) were also observed mostly at low abundance level. Some of the highly abundant star sequences observed by us were miR-106-b* (N1, TPM- 288.28); miR-17* (N1, TPM-1321.41); miR-92a-1* (N1, TPM- 942.78), miR- 25* (N1, TPM-265.44), and miR-374a* (N1, TPM-201.34). The results suggest that sRNA sequencing is a good approach for studying miRNA* abundance.

8miRNAs are often present in gene clusters and it is generally believed that individual miRNAs are generated from common polycistronic transcript by post - transcriptional processing. Since the genes in clusters are co-regulated and co-transcribed it is expected that the levels of these miRNAs are likely to be similar. This was investigated by checking the relative levels of individual miRNAs present in clusters using the expression patterns derived from deep sequencing. A total of 265 miRNA genes organized in 68 clusters were studied. As the sequencing reads represent mainly mature miRNAs and not individual precursors, it was not possible to analyze some of the clusters, such as those containing miR-15a, 16-1 and 15b, 16-2 that encode different precursors but the same mature miRNA. Furthermore, 19 clusters were also excluded from analysis as the miRNAs encoded in these clusters were either not expressed or expressed below 1 TPM based on our data. Therefore in this study results from analysis of only 20 clusters are shown (Additional file 2).

9The level of expression of different miRNAs in many clusters displayed variable expression as exemplified by the cluster containing miR-532 and 99b (Figure 5). In some cases the variation observed within a cluster was as much as 500 fold and the miRNA that showed highest level of expression was not the one closest to the transcription start site. For example, miR-25 showed highest expression though it was farthest from the start site in miR-106b cluster (Figure 5). This is likely to be due to variation in processing mechanisms at post - transcriptional level as has been seen in case of p53-mediated alteration in processing some of the miRNAs, such as miR-16-1. Overall the results suggest that the levels of miRNAs in any cell may be regulated by a number of different processes and mechanistic details of many of these are not yet known.

10Microarray based methods have been used extensively to profile expression levels of annotated genes at genomic scale. These methods cannot be used satisfactorily for analyzing expression levels of unknown transcripts. Deep sequencing of cellular RNAs is an alternative approach for deciphering expression profiles of genes. Though a few reports describing analysis of differentially expressed miRNAs deciphered by deep sequencing are available, there is no study yet on human PBMC. Moreover, there is no standardized tool available for profiling miRNA expression using deep sequencing of sRNA. Therefore we have studied normal peripheral blood mononuclear cells from two different individuals and cancer cells of myeloid lineage, K562 (chronic myelocytic leukemia) and HL60 (acute promyelocytic leukemia) by deep sequencing of sRNAs followed by analysis of the sequences using a custom designed computation pipeline. We used intersection of two methods for finding differentially expressed genes; a) SAM (significance analysis of microarrays) that allows for the control of false detection rate (FDR) and b) fold change, that is, miRNAs that showed more than 2.5 fold differences when compared to both normal samples. A list of differentially regulated miRNAs is shown (Figure 6). Interestingly, miR-1, 101, 106b, 146b-5p, 151-3p, 192, 21, 22, 27b, 30e and 361-3p displayed low expression in both the cell lines indicating a common role in leukemia genesis or progression (Figure 6, see inset). In K562 most of the miRNAs were down regulated except miR-486-3p and miR-504. The expression patterns derived from sRNA sequences were validated by RNase protection assay (RPA) which allowed both qualitative and quantitative analysis of RNA levels. Down regulation of expression of some of the miRNAs, such as miR-16, 22, 27a, 192 and let-7g in CML cell line K562 was observed using RPA assay (Figure 7). RNU6B was used as a loading control, for all the samples. We have also used quantitative RT-PCR for validating some of the results derived by deep sequencing and RPA. Fold differences in the expression of some of the differentially expressed genes in normal and K562 cells (miR-22 and miR-27a) using all three methods were found to be comparable indicating quantitative nature of these approaches (Additional file 3).

11While miR-21 was found to be one of the down regulated miRNAs in our study, it is expressed at high level in all the cancers tested till date indicating that miRNA expression is likely to be cell type and context dependent. Reduced expression of some of the miRNAs, such as miR-16, 151 and 142 was previously reported in CML. Up regulation of the polycistronic miR-17-92 cluster in CD34+ cells of CML patients was initially described by Venturini et al., 2007 and later was challenged by Agirre et al., 2008. According to the latter group there is no significant increase in the expression of this cluster in CML cells. We too did not observe any induction of this cluster in K562 cell line in agreement with the latter group.

12A significant number of differentially expressed miRNAs of HL60 (downregulated-miR-101, 126, 27b, 7; up regulated- let-7a, let-7d, miR-181a, -181a*, -181b and miR-199b) were mapped to chromosome 9 (Figure 8). Differential expression of miR-181a family (miR-181a, 181a* and 181b) was previously reported. Some of the differentially regulated miRNAs of K562, such as miR-101, miR-27b and miR-24 were also encoded by chromosome 9. The genomic locations of some of these miRNAs were distant from the ABL gene locus. Alterations in chromosome 9 are associated with a large number of diseases, particularly cancer. Therefore, there may be a link between alterations in chromosome 9 and differential expression of miRNAs.

13Identification of putative targets may help to understand the biological role of the differentially regulated miRNAs. In this study two different approaches were used for target identification. In the first approach a number of different software tools were employed for target prediction and in the second approach mRNA expression profiles were generated and the genes that showed inverse correlation with miRNA expression were identified. The second approach is based on a number of observations where target mRNA degradation by miRNAs were reported. For computational prediction stringent criteria, that is, only those genes that are predicted to be target of specific miRNAs by at least five of the eleven established target prediction programs (DIANA-microT, MicroInspector, miRanda, MirTarget2, miTarget, NBmiRTar, PicTar, PITA, RNA22, RNAhybrid, and TargetScan/TargertScanS) integrated by miRecords and also showed expression levels that are inversely correlated with the specific miRNA levels were considered. Putative targets of differentially regulated miRNAs are shown in Additional file 4 (K562) and Additional file 5 (HL60). Genes targeted by three or more miRNAs are shown in Table 1. Genes that are targets of multiple miRNAs are likely to be tightly regulated and may show graded response on the basis of expression of different miRNAs. MEIS2, up regulated in K562, was found to be the target of up to five different miRNAs, down regulated in both K562 and HL60. MEIS2 has been previously reported to be important for myeloid leukemogenesis. Similarly, SMAD7 a down regulated gene in HL60 is predicted to be a target of three miRNAs, up regulated in HL60. It is known that SMAD7 regulates SMAD and mitogen activated kinase (MAPKs) signaling and controls erythroid and megakaryocytic differentiation of erythroleukemia cells. TRIB2, multiply targeted by three up regulated miRNAs in HL60 has been reported to be present at low levels in AML.

14Many intron-encoded miRNAs are processed in DROSHA independent way (miRtrons). We have also analysed expression levels of differentially regulated intronic miRNAs. In order to check if differential expression of intronic miRNAs is due to transcriptional or post - transcriptional mechanisms, the levels of intronic miRNAs were compared with that of the host transcripts. Comparative analysis revealed positive correlations with respect to 8 out of 14 intronic miRNA - mRNA pairs (Table 2). Unlike miRNAs in clusters, the expression of some of the intronic miRNAs may be transcriptionally controlled.

15The analysis of microarray data revealed induction in the expression of some of the miRNA biogenesis genes (RNASEN, DGCR8, XPO5, RAN) in K562 cell line (Table 3). However, DICER1 was found to be down regulated (1.79 fold). This would result in a likely increase in the pre-miRNA population in cytosol but decrease in accumulation of mature miRNAs. This is consistent with the observation that K562 cells contain relatively less amount of miRNAs and that 23 out of 25 differentially expressed miRNAs are down regulated (see Figure 6). This observation was not due to artifact as two independent sequencing of sRNAs from K562 cells gave similar results (data not shown here). Low DICER and let-7 levels were also observed in lung adenocarcinoma. In contrast, a higher level of DICER1 in prostate adenocarcinoma accounted for up regulation of 39 of 45 differentially expressed miRNAs. In HL60 some of the components of biogenesis machinery (RNASEN and XPO5) were found to be up regulated.

16In principle deep sequencing of sRNAs should generate sequences from as yet unannotated regions of the genome. The analysis of the sequences through computational pipeline (Figure 9) showed a large number of unannotated sequences that are encoded by either intergenic or intronic regions. Since miRNAs are predominantly encoded by intergenic and intronic regions, these sequences were analyzed by a set of computational tools to identify putative novel miRNAs. Since predictions are based on identifying miRNA precursors, genomic regions (70 nucleotides) surrounding the sRNA sequences were extracted (see "Methods" for details). A total of 370 (357 major + 13 minor) novel miRNAs were predicted using computational pipeline (Additional file 6). The sequences and chromosomal locations of the predicted novel miRNAs from intronic and intergenic regions are listed in the supplementary table (Additional files 7 and 8). More than 95% of the novel miRNAs showed frequency counts lower than 10. Many novel miRNAs in our dataset showed presence of 1-11 isomiRs with varying frequencies. For example, frequencies ranging from 1-280 were observed for eight isomiRs of jnuhsa-204. Presence of isomiRs certainly strengthens classification of these molecules as novel miRNAs. These datasets were also analyzed using miRDeep, a published tool for miRNA analysis of deep sequencing data. The number of novel miRNAs predicted by miRDeep is much less than that using the pipeline described here (Additional files 7 and 8). This may be due to stringent criteria used by miRDeep that quite often misses known miRNAs. The expression of a few novel miRNAs, missed by miRDeep but predicted by "in-house tools", was validated through real-time PCR (for precursor miRNA) or RNase protection assay (for mature miRNAs) as shown in Figure 10 and 11. Interestingly miRNAs that were present as singleton (for example- jnuhsa-4-3p, jnuhsa-93-3p) were also detected using RNase protection assay (Figure 11). We also mapped predicted novel miRNAs to clusters based on the inter-miRNA distance of less than 20 kb (Figure 12). Some of the clusters contained only novel miRNAs while many others were part of previously known clusters. Based on miRDeep prediction a total of 12 novel miRNAs showed similarity in the seed region with the known miRNAs indicating that these miRNAs are likely to belong to the same family and thus may share common biology (Additional files 7 and 8).

17Our effort to understand the microRNAome of human leukocytes through deep sequencing technology has given us several interesting observations. These include discovery of novel miRNAs and a valuable list of differentially expressed miRNAs in chronic and acute myelogenous leukemia. The diversity in the miRNA wealth was realized through the detection of large number of novel miRNAs. Our major findings are a list of expressed miRNAs (534 known and 370 novel) in leukocytes and the discovery that the expression of miRNAs may be controlled by regulating post - transcriptional events, such as manipulating the level of DICER, an enzyme involved in biosynthesis of miRNAs in tumor cells (K562). Generally, singletons in the sequenced samples are not taken into consideration, as these are likely to occur due to potential sequencing error or transcription noise. Our experimental results show convincingly that such miRNAs are likely to be real and not a result of artifact of sequencing.

18The known miRNAs found to be most abundant in all the four samples were members of the let-7 family, 103, 185 and 320a, which were also reported in peripheral blood mononuclear cells using "Taqman microRNA assay". Therefore our results are in agreement with previous studies based on a different methodology. There are multiple mechanisms that are likely to regulate miRNA levels similar to that of mRNAs. These include both transcriptional and post - transcriptional processes. It is generally believed that levels of miRNAs are regulated transcriptionally. In this study an attempt was made to understand the role of different mechanisms in controlling miRNA levels by studying levels of individual miRNAs present in clusters and intronic miRNAs. It is expected that all miRNAs in a cluster belong to one transcription unit and are supposed to be synthesized as one long precursor. Since many miRNAs of a cluster showed variable expression patterns it is likely that the levels are controlled using post - transcriptional mechanisms. Therefore, regulation of post - transcriptional processing may be a preferred mechanism, instead of transcription control in case of miRNAs present in clusters. Studies describing the role of post - transcriptional mechanism in regulating levels of a number of miRNAs have been recently published. On the other hand our results regarding relative levels of intronic miRNA and the host mRNAs do indicate that the levels of many intronic miRNAs are regulated at transcriptional level.

19The isomiRs of different miRNAs were also frequently found suggesting that these do not arise due to rare events. As 5' and 3' ends of mature miRNAs are defined by processing from DROSHA and DICER, respectively, it is likely that isomiRs with extra nucleotides at the 5' end comes from variation in DROSHA processing while with extra bases at the 3' end comes from variation in DICER processing. Another possibility is processing of mature miRNA ends by yet an unidentified end-processing enzyme in the nucleus. Whether the extra nucleotides of isomiRs affect the target recognition needs to be evaluated.

20A number of differentially expressed miRNAs and their potential targets were identified by comparing miRNA and mRNA levels of normal PBMC versus that of two myeloid cancer cell lines K562 and HL60 along with computation approach. This strategy was used as degradation of many mRNAs was observed as a result of miRNA action. Both K562 and HL60 have been used as reference cell lines for studying the cellular and molecular events involved in the proliferation and differentiation of chronic myeloid or acute promyelocytic leukemia, respectively. Many miRNAs that were previously reported to be involved as tumor suppressors, such as let-7 g, miR-101, 16 and 192 were found to be down regulated in the leukemia cell lines. Many of the differentially regulated miRNAs, identified in this study have been reported to be involved in various cellular processes like cell cycle (miR-192), apoptosis (miR-16, -126, -98) differentiation (miR-27a, -181, -342, -223) DNA repair (miR-24, -210) metastasis, erythroid maturation (miR-22), erythropoiesis (miR-24) and hematopoiesis (miR-142, -181); highlighting their putative role in leukemogenesis or progression. Strikingly, majority of the differentially regulated miRNAs were found to be down regulated in K562 cell line. It is possible that low level of DICER, as revealed by gene expression profiling, is responsible for an overall reduction in the miRNA population. Global down regulation of miRNAs in mouse T cells using PCR and conventional sequencing approach was reported before. However, the cause of the down regulation was not clear, as mRNA profiling was not carried out. DICER has also been implicated in mouse T cell functions. Our studies strongly suggest that there may be a link between DICER and miRNA levels. Many of the genes, identified in this study as the potential targets of differentially regulated miRNAs are known to be involved in cancer through their effects on cell differentiation (CDK6, LIFR), apoptosis (PIM1) or hematopoiesis (GATA2, TAL1). Interestingly, majority of the targets of the cancer-associated miRNAs, such as -ARHGEF, CDK6, ETV5, GATA2, PIM1, LIFR, TAL1, PPARG, RANBP17, TFRC FOXP1 LPP, NCOA2 and NR4A3 have been previously shown to be associated with leukemogenesis.

21One of the features of this study is the computation pipeline that identifies novel miRNAs among sRNA sequences. This pipeline is different from other available pipelines as it has custom designed tools for identification of novel miRNAs that involve pre-processing of the sequences, an exhaustive elimination pipeline, folding and filtering to get novel predictions. The evidences supporting the predictions are identification of potential precursors by specific precursor identification tools, presence of star sequences, presence of isomiRs and in some cases experimental validation. The fact that even a single sequence read may not be a noise or error, has been validated by observing corresponding sRNA experimentally and indicates the power of our computation pipeline in identification of novel miRNAs. The total number of novel miRNAs identified from all the four samples is 370. Therefore, total number of miRNAs encoded by the human genome may be much larger than current estimate.

22In conclusion we list about 904 miRNAs that are expressed in normal and cancerous leukocytes, nearly 41% of these are novel. Our analysis shows that some of the novel miRNAs are likely to be clustered in the genome similar to many known clusters of miRNAs. Moreover, we have also identified a number of miRNAs that are differentially expressed in cancer cells studied by us. These are likely to be new markers/cause of carcinogenesis. In addition we show that miRNA levels can be regulated at post - transcriptional processing stage. Our results also show that singletons in deep sequencing reads are unlikely to be sequencing artifacts. In K562 cells low levels of miRNAs is correlated with reduced level of the enzyme DICER. In addition we have developed an improved automated computation pipeline for analysis of deep sequencing data to obtain quantitative profiles of known and novel miRNAs. Based on our observations, it is now possible to generate leukocyte specific miRNA arrays in order to study expression profiles of all miRNAs relevant to leukocytes. Similar approach can be used to generate other tissue specific miRNAome.

23The human chronic myeloid leukemia blast crisis cell line, K562 and acute promyelocytic leukemia cell line, HL60 were obtained from National Centre for Cell Sciences, Pune and maintained in RPMI 1640 and DMEM (Gibco) medium, respectively. The medium was supplemented with 10% FBS and penicillin-streptomycin and maintained at 37°C with 5% CO2 in incubator chamber. Buffy coat of healthy blood donors (N1 and N2) were collected from volunteers. Red cell lysis buffer (0.144M NH4Cl, 0.01M NH4HCO3) was added to buffy coat to lyse the remnant RBCs and pure WBC population was obtained by centrifugation at 3000 g.

24Total RNA isolation was carried out from peripheral blood and cell lines using TRIzol® Reagent (Invitrogen) as per manufacturer's instruction. RNA preparations were stored at - 80-°C till further use. Small RNA (sRNA) population was isolated by separating 10 μg of total RNA on denaturing polyacrylamide gel electrophoresis (PAGE) and cutting a portion of the gel corresponding to the size 18-30 nucleotides based standard oligonucleotide markers. Adapter (5') was ligated to sRNA population and ligated RNAs (40-60 nt) were purified by running on urea PAGE. This was followed by 3' adapter ligation and purification of adapter ligated RNAs (70-90 nt) in a similar manner. Modified sRNAs were reverse transcribed and then PCR amplified with adapter specific primers and the amplified cDNAs were finally purified on Urea PAGE to generate cDNA tag libraries for sequencing by illumina genome analyzer. The average number of sequencing reads was around 4.9 million.

25Four sRNA sequencing data comprising of peripheral blood leukocytes of two normal individuals (N1, N2) and tumor cells K562 and HL60 were obtained from Illumina fast track sequencing services. For each sample a Sequence file and a Tag file was provided. The Tag file comprises of unique sequences with their corresponding frequency. Tag Files are generated post alignment as a summary of Sequence Files and every Sequence File has a corresponding Tag File. Tag Files are generated to give the researcher an indication of most common to most rare sequences in the dataset. The numerical frequency of each sequence in the Sequence File for gene expression gives a true indication of relative expression of sequence transcripts. The unique sequences in the Tag files contain a 3' adaptor sequence (TCGTATGCCGTCTTCTGCTTG). The amount of the 3' adapter is variable and is dependent on the length of the sRNA. Each delivered sequence is 33/35 bases in length. A part of the adaptor sequence is seen in each sequence if the sRNA is shorter than 33/35 bases. This adapter and segments of it needs to be trimmed for proper alignment to the transcriptome/genome.

26i). Removal of the adaptor sequences: Since the sequence of the adaptor is known, a perl script was written to trim the adaptors.

27ii). Clustering and removal of redundancy after removal of the adaptor sequence: Although the tag file contained unique sequences, there were some, which after the removal of the adaptor were redundant. These identical sequences were represented once and their frequency was summed up.

28iii). Sequences less than 10 nucleotides in length were excluded

29iv). Conversion into fasta format: The final trimmed file was then converted into fasta format where the unique header (sequence identity) retained the information of the sequences length and frequency. The sequence ID comprised of a running number along with the length and frequency of that sequence.

30The sRNA sequences obtained were annotated against the known databases using the following protocol-

31List of the databases used for the annotation/elimination pipeline were:

32a) Mature miRNAs: from miRNA registry, release 14 (includes a total of 904 miRNAs [718 mature (major) +186 stars (minors)].

33b) ncRNAs: from Ensembl "Homo_sapiens.NCBI36.56.ncrna.fa" (includes the precursor miRNAs and other ncRNAs like sn/sno/sca RNAs, tRNAs, rRNAs).

34c) RNA database from the FTP site NCBI (includes rRNAs and mRNAs).

35d) The exons were obtained from the Contig files by a self written script.

36e) Intergenic/intronic sequences: obtained through in-house built Perl script (using the Homo Sapiens Contig file 29 Feb, 2008 version). These sequences served as a source for finding novel miRNAs (intergenic/intronic).

37One of our objectives was to study the expression pattern of the known miRNAs. To generate the expression profile of the known miRNAs, the sRNA sequences of all the 4 samples were matched against the known miRNA sequences using BLASTN and an in house built pattern matching tool. The parameters used for BLAST were tuned to obtain maximum matches, such as the word size was set to 7 nucleotides, filtering was turned off and the number of alignments reported were increased.

38An in house - built script was written to do a fast matching of the sequences with the created databases. A mismatch of up to 2 nucleotides was allowed. The pool of unmatched sequences at the end of the pipeline served as a source of novel miRNAs (Figure 1).

39Normalisation was carried out as the total number of reads from different experiments was not same and variations in the number of reads of individual miRNA can be due to sequencing depth. The number of reads of a transcript/sequence (representing a known miRNA) was divided by the total clone count of the sample and multiplied by 106. The total clone count is the sum of the frequencies of all the unique sequences/transcripts present in the trimmed file. The Additional file 1 contains a list of the known miRNAs found in all the 4 samples along with their respective frequencies and TPM values.

40Differentially expressed miRNAs were identified by using a combination of two methods:

41We performed the t-test procedure within significance analysis of microarrays (SAM) to look for differentially expressed miRNAs. SAM calculates a score for each gene on the basis of the change in expression relative to the standard deviation of all measurements.

42The known miRNAs in normal samples were compared to the cancer cell line. The miRNAs showing more than 2.5 fold difference as compared to both the normal cells were considered as differentially regulated. The samples with less than 10 TPM in both normal and cancer samples were ignored.

43The strategy is based on first removing all known RNAs including those derived from exonic regions and then identifying those that are derived from intronic and intergenic regions. These were then subjected to some of the ab initio miRNA prediction algorithms like SCFG based CID-miRNA and CSHMM that use stringent criteria to distinguish between real and pseudo miRNAs. Besides this we ensure high discriminative power by using the following filters:

44a) Occurrence of the sequence in the stem region of the precursor

45b) Presence of IsomiRs

46c) Presence of minor/star sequences.

47d) Taking a consensus among the prediction tools

48Following is the detailed pipeline used for the prediction of novel miRNAs:

49The unmatched sequences (from the elimination pipeline) were matched to the intergenic/intronic regions. The exact matched sequences were extracted along with 70 nucleotides flanking both the ends representing potential precursor sequences.

50The sequences were scanned for presence of potential precursor miRNA using CID-miRNA and CSHMM prediction tools. The folded sequences generated were then checked to see if the sRNA (putative mature miRNA obtained by sequencing) occurs in the folded putative precursor as the window scanning approach used could report a folded structure not involving the concerned sRNA. Only those hairpins were kept which contained sRNA. The next step involved locating the position of the sRNA in the hairpins. Since mature miRNAs are known to be arising from the stem portion and not the loop, only those hairpins in which the sRNAs occurred in the stem were classified as correct cases and the remaining as prediction errors. These correct cases were further tested by MiPred.

51A list of all the predicted correct precursor sequences was created and the sRNAs derived from common precursors were grouped into a common family. The sRNAs derived from the same precursors were kept together in a family. The most abundant member was designated as the mature miRNA. The sRNAs that differed from the representative by a few nucleotides were called IsomiRs and those that had a different, partially complementary sequence and were located in the other strand (stem of the hairpin loop) were called stars. The Additional files 7 and 8 comprise of the novel miRNAs grouped into families on the basis of sRNAs falling within the same precursor. The representative sequence was chosen on the basis of abundance. The most abundant sequence was selected as representative of the family. The scores from the 4 tools (CID-miRNA, CSHMM, miRDeep, MiPred assigned to the corresponding precursors are also listed).

52On the basis of the positions of the predicted novel miRNAs, their tendency to occur in a known cluster or a new cluster was checked. miRNAs located within 20 kb of the known or novel miRNA were considered as part of the same cluster.

53The novel miRNAs from the 4 samples (both intergenic and intronic) were pooled and the redundancy was removed to get a final set of novel miRNAs. These were given a unique name (Additional file 6). The representative miRNAs in the Additional files 7 and 8 along with their sample IDs also have these names.

54The RPA assay for mature miRNA detection was done using mirVana ™ miRNA detection kit as per supplier's instructions (Applied Biosystems). To detect the precursor miRNA, total RNA was treated with DNAase I (MBI fermentas). Reverse transcription was done using superscript® III (Invitrogen) reverse transcriptase as per supplier's instructions. A list of primers used for real time RT-PCR amplification of precursors of novel miRNAs is given in Additional file 9. Real-time PCR was done with following parameters- Initial denaturation - 94°c for 2 min, denaturation- 94°c for 30 sec, annealing- 60°c for 1 min for 40 cycles using SYBR Green PCR Master Mix (Applied biosystems).

55Total RNA from normal PBMC and two myeloid leukemia cancer samples (K562 and HL60) was sent for gene expression profiling using Ocimum microarray platform (Ocimum, Hyderabad, India). The expression data for each sample was generated on Affymetrix Human Gene 1.0 ST arrays. A fold change cut-off of ± 1.5 resulted into FDR < 0.001 for both the comparisons. The quality control analysis was carried out using Affymetrix Expression Console (EC). The statistical analysis was performed using R-programming language and the biological analysis was carried out using GenowizTM software. The data obtained was normalized and genes showing more than 1.5 fold differences in the cancer cell lines as compared to normal were marked as differentially expressed genes. A total of 1856 genes were found to be up regulated and 1696 were down regulated in K562. In HL60, 1497 genes were up regulated and 1213 genes were found as downregulated (Additional file 10).

56The most probable targets of the differentially regulated miRNAs were fished out using following two criteria - 1. Prediction by at least five of the established target prediction programs- A list of putative targets was prepared (List A) using intersection between at least five of the eleven established target prediction programs compiled by mIRecords: (DIANA-microT, MicroInspector, miRanda, MirTarget2, miTarget, NBmiRTar, PicTar, PITA, RNA22, RNAhybrid, and TargetScan/TargertScanS). 2. Inverse correlation in expression pattern between miRNA and coding genes- The putative target genes in list A was compared to the list of differentially regulated genes (showing more than 1.5 fold difference) (List B) in K562 and HL60 as obtained from expression profiling data, and only those genes that show inverse correlation to the miRNA levels were considered as most genuine putative targets of the select miRNAs.

57The entire computational pipeline described in this paper is available at the website: http://mirna.jnu.ac.in/deep_sequencing/deep_sequencing.html. The software is also available on request.

58The data discussed in this publication have been deposited in NCBI's Gene Expression Omnibus and are accessible through GEO Series accession number GSE19833 Super series: http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE19833

59The GSE19833 Super series comprises of the following sub series:

60Microarray Expression data: GSE 19789 (consists of raw data and normalized data) Deep sequencing data: GSE 19812 sequence data (consists of the raw sequences, untrimmed tag and trimmed tag files and the expression profile of the known miRNAs (Additional File 1).

61miRNA: microRNA; UTR: untranslated region; RISC: RNA induced silencing complex; PTEN: phosphatase and tensin homolog; RAS: rat sarcoma viral oncogene homolog; HMGA2: high mobility group AT-hook 2; BCL2: B-cell CLL/lymphoma 2; PBMC: peripheral blood mononuclear cell; SAM: significance analysis of microarrays; FDR: false discovery rate; CML: chronic myelogenous leukemia; AML: myelogenous leukemia; ABL: c-abl oncogene 1, receptor tyrosine kinase; MEIS2: Meis homeobox 2; SMAD7: SMAD family member 7; TRIB2: tribbles homolog 2 (Drosophila); MAPK: mitogen activated protein kinase; BCL11A: B-cell CLL/lymphoma 11A (zinc finger protein); ETV4: ets variant 4; IDH1: isocitrate dehydrogenase 1 (NADP+), soluble; PPARG: peroxisome proliferator-activated receptor gamma; EIF4A2: eukaryotic translation initiation factor 4A, isoform 2; LCK: lymphocyte-specific protein tyrosine kinase; NCOA2: nuclear receptor coactivator 2; NR4A3: nuclear receptor subfamily 4, group A, member 3; PIM1: proviral integration site 1; RNASEN: ribonuclease type III, nuclear; DGCR8: DiGeorge syndrome critical region gene 8; XPO5: exportin 5; RAN: ras-related nuclear protein; PRMT5: protein arginine methyltransferase 5; FXR2: fragile X mental retardation, autosomal homolog 2; MOV10: Moloney leukemia virus 10; SNP: single nucleotide polymorphism

62AB and LK conceptualized the study. CV generated the computational pipeline, processed, normalized and annotated the data and did the novel miRNA prediction analysis. RK and HMA conducted the miRNA differential expression, miRNA cluster, microarray data and target prediction analyses and wet experiments. PS provided the RNA samples for sequencing. RG conducted SAM analyses. RK and AB wrote the manuscript. All the authors read and approved the final manuscript.