Step 1: Process of H3K4me3 and Pol II data
peak_merge(): Merge one kind of peaks
(H3K4me3 or Pol II)
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.
library(primirTSS)
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
peak_join(): Join two kinds of peaks
(H3K4me3 and Pol II)
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
find_tss is the primary function in the package. The
program will first score the candidate TSSs of miRNA and pick up the
best candidate in the first step of prediction, (where users can set
flanking_num and threshold).- After the first step, H3K4me3 and Pol II data, miRNA expression
profiles and DHS check, protein-coding genes expression profiles, if
provided, will be integrated to decrease the rate of false
positive.
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.
ignore_DHS_check: If users do not have their own miRNA
expression profile, the function will employ all the miRNAs already
annotated in human, but we suggest using DHS data of the cell line from
ENCODE to check whether this miRNA is expressed in the
cell line or not as well as and all human gene expression profiles from
Ensemble to check the relative position of TSSs and
protein-coding genes to improve the accuracy of prediction.
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.
expressed_mir: If users have their own miRNA expression
profiles, we will use the expressed miRNAs and we suggest not using DHS
data of the cell line or others to check the expression of miRNAs.But
the protein-coding gene check to check the relative position of TSSs and
protein-coding genes is necessary, which helps to verify the precision
of prediction.
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)
expressed_gene: Additionally, users can also specify
certain genes expressed in the cell-line being analyzed:
Step4: Analysis of results
Here is a demo of predicting TSS for hsa-mir-5697, ignore DHS
check.
PART1, $tss_df:
ownmiRNA$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:
| 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:
host_TSS: The TSSs of miRNA that are close to
the TSS of protein-coding gene implying they may share the same TSS, on
the condition where mir_context = intra. (See above:
mir_context)
intra_TSS: The TSSs of miRNA that are NOT close
to the TSS of protein-coding gene, on the condition where
mir_context = intra.
overlap_inter_TSS: The TSSs of miRNA are
cataloged as overlap_inter_TSS when the pri-miRNA gene
overlaps with Ensembl gene, on the condition where
“mir_context = inter”.
inter_inter_TSS: The TSSs of miRNA are cataloged
as inter_inter_TSS when the miRNA gene does NOT overlap
with Ensembl gene, on the condition where
“mir_context = inter”.
(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:
find_nearest_peak_log: If no peaks
locate in the upstream of a stem-loop to help determine putative TSSs of
miRNA, we will fail to find the nearest peak, and this miRNA will be
logged in find_nearest_peak_log.
eponine_score_log: For a certain
miRNA, if none of the candidate TSSs scored with Eponine method meet the
threshold we set, we will fail to get an eponine score, and this miRNA
will be logged in eponine_score_log.
DHS_check_log: For a certain miRNA,
if no DHS signals locate within 1 kb upstream of each putative TSSs,
these putative TSSs will be filtered out, and this miRNA will be logged
in DHS_check_log.
gene_filter_log: For a certain
miRNA, when integrating expressed_gene data to improve prediction, if no
putative TSSs are confirmed after considering the relative position
relationship among TSSs, stem-loops and expressed genes, this miRNA will
be filtered out and logged in gene_filter_log.
Figure S1. Visualized result for miRNA TSSs by
Figure S2. Graphical web interface of
Figure S3. Graphical web interface of