1 Introduction

Identifying human miRNA transcriptional start sites (TSSs) plays a significant role in understanding the transcriptional regulation of miRNA. However, due to the quick capping of pri-miRNA and many miRNA genes may lie in the introns or even exons of other genes, it is difficult to detect miRNA TSSs. miRNA TSSs are cell-specific. And miRNA TSSs are cell-specific, which implies the same miRNA in different cell-lines may start transcribing at different TSSs.

High throughput sequencing, like ChIP-seq, has gradually become an essential and versatile approach for us to identify and understand genomes and their transcriptional processes. By integrating H3k4me4 and Pol II data, parting of false positive counts after scoring can be filtered out. Besides, DNase I hypersensitive sites(DHS) also imply TSSs, where miRNAs will be accessible and functionally related to transcription activities. And additionally, the expression profile of miRNA and genes in certain cell-line will be considered as well to improve fidelity. By employing all these different kinds of data, here we have developed the primirTSS package to assist users to identify miRNA TSSs in human and to provide them with related information about where miRNA genes lie in the genome, with both command-line and graphical interfaces.

2 Find the best putative TSS


1 primirTSS

Install the latest release of R, then get primirTSS by starting R and entering the commands:

if (!requireNamespace("BiocManager", quietly = TRUE))

Or install the development version of the package using the same command, but with GitHub account name.


2 Install Java SE Development Kit(JDK)

As Java development environment is indispensable for the primary function in our package, it is necessary for users to install Java SE Development Kit 10 before using primirTSS.

3 Load the package into R session

#> Warning: replacing previous import 'utils::findMatches' by
#> 'S4Vectors::findMatches' when loading 'phastCons100way.UCSC.hg38'

3 Getting Started

Step 1: Process of H3K4me3 and Pol II data

H3K4me3 and Pol II data are key points for accurate prediction with our method. If one of these two peak data is input, before executing the main function find_TSS, the function peak_merge should be used to merge adjacent peaks whose distance between each other is less than n base pairs and return the merged peaks as an output.

peak_df <- data.frame(chrom = c("chr1", "chr2", "chr1"),
                       chromStart = c(450, 460, 680),
                       chromEnd = c(470, 480, 710),
                       stringsAsFactors = FALSE)
peak <-  as(peak_df, "GRanges")
peak_merge(peak, n =250)
#> GRanges object with 2 ranges and 0 metadata columns:
#>       seqnames    ranges strand
#>          <Rle> <IRanges>  <Rle>
#>   [1]     chr1   450-710      *
#>   [2]     chr2   460-480      *
#>   -------
#>   seqinfo: 2 sequences from an unspecified genome; no seqlengths

If both of H3K4me3 and Pol II data, after separately merging these two kinds of evidence first, peak_join should be employed to integrate H3K4me3 and Pol II peaks and return the result as bed_merged parameter for the main function find_tss.

peak_df1 <- data.frame(chrom = c("chr1", "chr1", "chr1", "chr2"),
                       start = c(100, 460, 600, 70),
                       end = c(200, 500, 630, 100),
                       stringsAsFactors = FALSE)
peak1 <-  as(peak_df1, "GRanges")

peak_df2 <- data.frame(chrom = c("chr1", "chr1", "chr1", "chr2"),
                       start = c(160, 470, 640, 71),
                       end = c(210, 480, 700, 90),
                       stringsAsFactors = FALSE)
peak2 <-  as(peak_df2, "GRanges")

peak_join(peak1, peak2)
#> GRanges object with 3 ranges and 0 metadata columns:
#>       seqnames    ranges strand
#>          <Rle> <IRanges>  <Rle>
#>   [1]     chr1   160-200      *
#>   [2]     chr1   470-480      *
#>   [3]     chr2     71-90      *
#>   -------
#>   seqinfo: 2 sequences from an unspecified genome; no seqlengths

Step 2: Predict most possible TSS for miRNA

There will be different circumstances where not all miRNA expression profiles, DHS data, protein-coding gene(‘gene’) expression profiles are available:

Circumstance 1: no miRNA expression data; then suggest DHS check and protein-coding gene check.

peakfile <- system.file("testdata", "HMEC_h3.csv", package = "primirTSS")
DHSfile <- system.file("testdata", "HMEC_DHS.csv", package = "primirTSS")
peak_h3 <- read.csv(peakfile, stringsAsFactors = FALSE)
DHS <- read.csv(DHSfile, stringsAsFactors = FALSE)
DHS <- as(DHS, "GRanges")
peak_h3 <-  as(peak_h3, "GRanges")
peak <- peak_merge(peak_h3)
no_ownmiRNA <- find_tss(peak, ignore_DHS_check = FALSE,
                        DHS = DHS, allmirdhs_byforce = FALSE,
                        expressed_gene = "all",
                        allmirgene_byforce = FALSE,
                        seek_tf = FALSE)

Circumstance 2: miRNA expression data provided; then no need for DHS check but protein-coding gene check.

bed_merged <- data.frame(
                chrom = c("chr1", "chr1", "chr1", "chr1", "chr2"),
                start = c(9910686, 9942202, 9996940, 10032962, 9830615),
                end = c(9911113, 9944469, 9998065, 10035458, 9917994),
                stringsAsFactors = FALSE)
bed_merged <- as(bed_merged, "GRanges")

expressed_mir <- c("hsa-mir-5697")

ownmiRNA <- find_tss(bed_merged, expressed_mir = expressed_mir,
                     ignore_DHS_check = TRUE,
                     expressed_gene = "all",
                     allmirgene_byforce = TRUE,
                     seek_tf = FALSE)

Step 3: Searching for TFs

Step4: Analysis of results

Here is a demo of predicting TSS for hsa-mir-5697, ignore DHS check.

PART1, $tss_df:

#> # A tibble: 1 × 10
#>   mir_name     chrom stem_loop_p1 stem_loop_p2 strand mir_context tss_type gene 
#>   <chr>        <chr>        <dbl>        <dbl> <chr>  <chr>       <chr>    <chr>
#> 1 hsa-mir-5697 chr1       9967381      9967458 +      intra       host_TSS ENSG…
#> # ℹ 2 more variables: predicted_tss <dbl>, pri_tss_distance <dbl>

The first part of the result returns details of predicted TSSs, composed of seven columns: mir_name, chrom, stem_loop_p1, stem_loop_p2, strand mir_context, tss_type gene and predicted_tss:

Entry Implication
mir_name Name of miRNA.
chrom Chromosome.
stem_loop_p1 The start site of a stem-loop.
stem_loop_p2 The end site of a stem-loop.
strand Polynucleotide strands. (+/-)
mir_context 2 types of relative position relationship between stem-loop and protein-coding gene. (intra/inter)
tss_type 4 types of predicted TSSs. See the section below TSS types for details.(host_TSS/intra_TSS/overlap_inter_TSS/inter_TSS)
gene Ensembl gene ID.
predicted_tss Predicted transcription start sites(TSSs).
pri_tss_distance The distance between a predicted TSS and the start site of the stem-loop.

TSSs are cataloged into 4 types as below:

(See Xu HUA et al 2016 for more details)

PART2, $log:

The second part of the result returns 4 logs created during the process of prediction:

4 Plot the prediction of TSS for miRNA

plot_primiRNA(expressed_mir, bed_merged,
              flanking_num = 1000, threshold = 0.7,
              ignore_DHS_check = TRUE,
              DHS, allmirdhs_byforce = TRUE,
              expressed_gene = "all",
              allmirgene_byforce = TRUE)

plot_tss    Figure S1. Visualized result for miRNA TSSs by Plot pri-miRNA TSS()

As Figure S1 shows, the picture contains information of the pri-miRNA’s coordinate, the closest gene to the miRNA, the eponine score of the miRNA’s candidate TSS and the conservation score of the miRNA’s candidate TSS. There are six tracks plotted in return:

Entry Implication
Chromosome Position of miRNA on the chromosome.
hg19 Reference genome coordinate in hg19.
pri-miRNA: Position of pri-miRNA.
Ensemble genes Position of related protein-coding gene.
Eponine score Score of best putative TSS by Eponine method.
Conservation score Conservation score of TSS.                         

5 Graphical web interface for prediction

TAG1: Find the best putative TSS

shiny_tss  Figure S2. Graphical web interface of Find pri-miRNA TSS()

As Figure S2 shows, if we want to use the shiny app, we should select the appropriate options or upload the appropriate files. Histone peaks, Pol II peaks and DHS files are comma-separated values (CSV) files, whose first line is chrom,start,end. Every line of miRNA expression profiles has only one miRNA name which start with hsa-mir, such as hsa-mir-5697. Every line of gene expression profiles has only one gene name which derived from Ensembl, such as ENSG00000261657. All of miRNA expression profiles and gene expression profiles do not have column names. If we have prepared, we can push the Start the analysis button to start finding the TSSs. The process of analysis may need to take a few minutes, and a process bar will appear in right corner.

As a result, we will view first six rows of the result. The first five columns are about miRNA information, next five columns are about TSS information. The column of gene denotes the gene whose TSS is closest to the miRNA TSS. The column of pri_tss_distance denotes the distance between miRNA TSS and stem-loop. If users choose to get TFs simultaneously, they will have an additional column, tf, which stores related TFs.

TAG2: Plot pri-miRNA TSS

shiny_plot  Figure S3. Graphical web interface of Plot pri-miRNA TSS()

As Figure S4 shows, if we select the appropriate options and upload the appropriate files, we can have a picture of miRNA TSSs.

Session info

Here is the output of sessionInfo() on the system on which this document was compiled:

#> R version 4.4.0 beta (2024-04-15 r86425)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 22.04.4 LTS
#> Matrix products: default
#> BLAS:   /home/biocbuild/bbs-3.19-bioc/R/lib/ 
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> time zone: America/New_York
#> tzcode source: system (glibc)
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> other attached packages:
#> [1] primirTSS_1.22.0
#> loaded via a namespace (and not attached):
#>   [1] tidyselect_1.2.1                  dplyr_1.1.4                      
#>   [3] blob_1.2.4                        filelock_1.0.3                   
#>   [5] R.utils_2.12.3                    Biostrings_2.72.0                
#>   [7] bitops_1.0-7                      fastmap_1.1.1                    
#>   [9] RCurl_1.98-1.14                   BiocFileCache_2.12.0             
#>  [11] GenomicAlignments_1.40.0          XML_3.99-0.16.1                  
#>  [13] digest_0.6.35                     lifecycle_1.0.4                  
#>  [15] KEGGREST_1.44.0                   RSQLite_2.3.6                    
#>  [17] magrittr_2.0.3                    compiler_4.4.0                   
#>  [19] rlang_1.1.3                       sass_0.4.9                       
#>  [21] tools_4.4.0                       utf8_1.2.4                       
#>  [23] yaml_2.3.8                        BSgenome.Hsapiens.UCSC.hg38_1.4.5
#>  [25] rtracklayer_1.64.0                knitr_1.46                       
#>  [27] S4Arrays_1.4.0                    phastCons100way.UCSC.hg38_3.7.1  
#>  [29] bit_4.0.5                         curl_5.2.1                       
#>  [31] DelayedArray_0.30.0               abind_1.4-5                      
#>  [33] BiocParallel_1.38.0               HDF5Array_1.32.0                 
#>  [35] withr_3.0.0                       purrr_1.0.2                      
#>  [37] R.oo_1.26.0                       BiocGenerics_0.50.0              
#>  [39] grid_4.4.0                        stats4_4.4.0                     
#>  [41] fansi_1.0.6                       Rhdf5lib_1.26.0                  
#>  [43] SummarizedExperiment_1.34.0       cli_3.6.2                        
#>  [45] rmarkdown_2.26                    crayon_1.5.2                     
#>  [47] GenomicScores_2.16.0              generics_0.1.3                   
#>  [49] httr_1.4.7                        rjson_0.2.21                     
#>  [51] DBI_1.2.2                         cachem_1.0.8                     
#>  [53] rhdf5_2.48.0                      stringr_1.5.1                    
#>  [55] zlibbioc_1.50.0                   parallel_4.4.0                   
#>  [57] AnnotationDbi_1.66.0              BiocManager_1.30.22              
#>  [59] XVector_0.44.0                    restfulr_0.0.15                  
#>  [61] matrixStats_1.3.0                 vctrs_0.6.5                      
#>  [63] Matrix_1.7-0                      jsonlite_1.8.8                   
#>  [65] IRanges_2.38.0                    S4Vectors_0.42.0                 
#>  [67] bit64_4.0.5                       tidyr_1.3.1                      
#>  [69] jquerylib_0.1.4                   glue_1.7.0                       
#>  [71] codetools_0.2-20                  JASPAR2018_1.1.1                 
#>  [73] stringi_1.8.3                     BiocVersion_3.19.1               
#>  [75] GenomeInfoDb_1.40.0               BiocIO_1.14.0                    
#>  [77] GenomicRanges_1.56.0              UCSC.utils_1.0.0                 
#>  [79] tibble_3.2.1                      pillar_1.9.0                     
#>  [81] rappdirs_0.3.3                    htmltools_0.5.8.1                
#>  [83] rhdf5filters_1.16.0               GenomeInfoDbData_1.2.12          
#>  [85] BSgenome_1.72.0                   R6_2.5.1                         
#>  [87] dbplyr_2.5.0                      evaluate_0.23                    
#>  [89] lattice_0.22-6                    Biobase_2.64.0                   
#>  [91] AnnotationHub_3.12.0              png_0.1-8                        
#>  [93] R.methodsS3_1.8.2                 Rsamtools_2.20.0                 
#>  [95] memoise_2.0.1                     bslib_0.7.0                      
#>  [97] SparseArray_1.4.0                 xfun_0.43                        
#>  [99] MatrixGenerics_1.16.0             pkgconfig_2.0.3