yue-wang-biomath/RgnTX
Colocalization analysis of transcriptome elements in the presence of isoform heterogeneity and ambiguity

Package index
Search the yue-wang-biomath/RgnTX package
Vignettes
Functions
63
Source code
35
Man pages
31
Home
/
GitHub
/
yue-wang-biomath/RgnTX
/

In yue-wang-biomath/RgnTX: Colocalization analysis of transcriptome elements in the presence of isoform heterogeneity and ambiguity 

library(knitr)
opts_chunk$set(concordance=FALSE)
set.seed(12345677)

library(RgnTX)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)

Introduction

Recent development of high-throughput sequencing technology has generated large amounts of data at transcriptome level such as RNA-binding protein target sites, RNA modifications and other RNA-related genomic features. These transcriptome attributes can be converted to a set of locations or regions.

Currently, colocalization analysis of two genomic region sets has been widely used to infer potential interactions between corresponding biological attributes. A number of methods have been developed for this. For example, r BiocStyle::Biocpkg("regioneR") has been created to statistically estimate association between genomic region sets using a permutation test framework. All its functions are genome and mask-aware. However, there has been no such tools enable researchers to make data analysis and interpretation at transcriptome level. Here RgnTX is invented to fill this blank.

Inspired by r BiocStyle::Biocpkg("regioneR"),RgnTX extends an assumption-free permutation test framework transcriptome-wide. Different from existing approaches, RgnTX allows the integration of transcriptome annotations so as to model the complex alternative splicing patterns. Importantly, it also supports the testing of transcriptome elements without clear isoform association, which is often the real scenario due to technical limitations.

RgnTX makes a useful tool for transcriptome region set association analysis. In this vignette, with some simulated datasets, we show that permutation tests conducted in the genome and the transcriptome are significantly different and may return distinct results in section 2. The most core functions for conducting permutation tests and examples with real datasets are introduced in section 6.

What is the difference between permutation tests in the genome and the transcriptome?

What does a permutation test do?

Suppose we want to know if a region set A has some association with another region set B. We could pick some random region sets, and calculate the number of overlaps between them and B. If the overlap numbers of these random region sets and B are generally smaller than that of A and B, it is more likely to say there exists some association between A and B, otherwise not.

Such permutation test framework mainly involves two steps: the randomization step (picking random region sets) and the evaluation step (evaluating overlapping counts or other test statistic). Importantly, performing it over the transcriptome and the genome can be different in both of these two steps.

Randomization step

Compared with the linear genome space, transcriptome space is heterogeneous due to the existence of multiple isoform transcripts, the splicing junctions and exons used by different frequencies. The isoform specificity of transcriptome elements is often lost in the process of mapping due to technical limitations. To perform isoform-aware permutation for transcriptome element with isoform ambiguity, RgnTX considers only the transcripts that overlap with it when projected to the genome, so as to retrain maximal isoform information. In this example shown by the following figures, feature 3 will be permutated on transcript Tx3 only, while feature 2 will be permutated on all the 3 transcripts.

knitr::include_graphics("figures/section2.1.png")

knitr::include_graphics("figures/section2.2.png")

knitr::include_graphics("figures/section2.3.png")

This is supported by the RgnTX function randomizeFeaturesTx and the function randomizeFeaturesTxIA. The former is for features with known isoform belongings. The latter is for features with isoform ambiguity (IA). Besides, RgnTX also provides function randomizeTx supporting picking random region sets within any specified transcriptome space. Users can choose a suitable strategy to do randomization conveniently within only a few lines of codes. More details about this are introduced in section 5.

Evaluation step

Another non-trivial task is to develop evaluation function for a test statistic, such as the number of overlaps, to be tested between two region sets. A key limitation of existing evaluation approaches is that, they were designed primarily for genome-based evaluation. Two transnscriptome elements with shared exons but are located on the independent transcripts are also considered to have overlaps by them. The information of transcripts is totally lost during comparison, which can induce substantial bias when applying to analyze transcriptome region sets.

RgnTX provides several strategies for transcriptome-level evaluation. Function overlapCountsTX counts the number of overlaps between two transcriptome elements in the transcriptome directly. Function overlapCountsTXIA calculates a weighted number of overlaps between a feature set (with isoform ambiguity) and a transcriptome region set (without isoform ambiguity). It is also possible to write a custom evaluation function and pass it to the main permutation test function. Details about this part are introduced in section 7.

Comparison of permutations over different kinds of spaces

We perform a series of simulations to show the difference of permutation tests conducted in the transcriptome and in the genome.

In this simulation, we randomly generated 1000 pairs of region sets in the following four spaces, with each region set containing 500 regions of 200nt long on the same 300 genes, and then count the number of overlaps between each pair.

Figure \@ref(fig:Fig4) shows the evaluation results. In the first four cases, overlaps are counted based on the genome (or after projection). In the last case, overlaps are counted in the transcriptome directly. The former four cases are evaluated by r BiocStyle::Biocpkg("regioneR") function numOverlaps, while the last is evaluated by RgnTX function overlapCountsTx.

knitr::include_graphics("figures/section2.4.png")

As the results show, permutations over different genome and transcriptome space return distinct results. A small number of overlaps are observed in the transcriptome (the last case), which is not surprising due to the existence of multiple isoform transcripts. Interestingly, the results from exonic DNA and mRNA (projected to genome) are also different, which is due to exons used at different frequencies by isoform transcripts. These results suggest that the permutation of heterogeneous transcriptome elements can be substantially more complex than genome-based elements. Genome-based methods are not appropriate to analyze transcriptome elements. That’s why we develop RgnTX to do this.

Installation

To install this package via devtools, please use the following codes.

if (!requireNamespace("devtools", quietly = TRUE))
    install.packages("devtools")

devtools::install_github("yue-wang-biomath/RgnTX", build_vignettes = TRUE)

To install this package via BiocManager, please use the following codes.

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("RgnTX")

Feature/region format in RgnTX

RgnTX is based on the r BiocStyle::Biocpkg("GenomicRanges") and r BiocStyle::Biocpkg("GenomicFeatures") structure. In this section, we introduce the format of region sets and other associated metadata used in RgnTX.

A TxDb object stores transcripts, exons, CDS and other types of genomic features. A TxDb object is usually required by RgnTX functions. Currently, all the examples in this vignette are based on the TxDb.Hsapiens.UCSC.hg19.knownGene. Users can assign other TxDb objects via related function arguments.

A TxDb object provides default ids and names for transcripts. According to r BiocStyle::Biocpkg("GenomicFeatures"), a transcript is guaranteed to be related to a unique transcript id. But a transcript's name is not guaranteed to be unique or even defined. In this vignette, all the examples use transcript ids instead of names to indicate specific transcripts.

Please pay attention not to confuse the transcript ids with other types of ids, such as ids for CDS, exons and other genomic features provided by a TxDb object, or external ids for genes like Entrez Gene and Ensembl ids. Functions in RgnTX will not work with wrong types of ids.

The feature/region sets in RgnTX should be either GRanges or GRangesList.

For a GRangesList feature set, every list element that contains a set of ranges is a feature. For a GRanges feature set without special instructions, a single range is a feature. For a GRanges feature set containing a metadata column named "group", the set of ranges having the same group number is a feature.

The transcript ids (if available) should be provided as follows. For GRangesList objects, the name of each list element (feature) should be the id of the transcript that this feature is located on. For GRanges objects, transcript ids should be contained in a metadata column named "transcriptsHits".

Here is an example of GRangesList feature set following the format required by RgnTX. Its names should be transcript ids.

library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene

set.seed(12345677)
example.list <- randomizeTx(txdb, random_num = 10, random_length = 200)
example.list

The following GRanges object represents the same region set as the above one. Ranges having the same group number represent a feature. The transcriptsHits column contains corresponding transcript ids. Both formats are accepted by RgnTX functions and can be transformed to each other via function GRangesList2GRanges and function GRanges2GRangesList.

example.gr <- GRangesList2GRanges(example.list)
example.gr
example.list <- GRanges2GRangesList(example.gr)

Randomization function

In RgnTX, there are mainly two kinds of randomization functions. One can pick random regions over a specified space (see section 5.1). One can also randomize specified features into transcriptome, i.e., input a set of features and get a randomized result of it (see section 5.2).

Pick random regions over a specified space

In this strategy, users need to specify a scope that random region sets will be picked from.

randomizeTx function

In this function, users can pass transcript ids via an argument to make random regions to be picked from the transcripts having these ids.

Step 1: specify a randomization space

One needs to determine three arguments: a TxDb data, a set of transcript ids and the type of randomized regions the function will return. As an example, we specify a randomization space formed by five mRNAs.

txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene 
type <- "mature"
trans.ids <- c("170", "782", "974", "1364", "1387")

We feed them into the randomizeTX function.

Step 2: randomizeTx function

We use this function to pick 10 random regions, and each region is 100 nt long. Each region is picked from a random transcript among the input transcripts.

set.seed(12345677)
randomResults <- randomizeTx(txdb, trans_ids = trans.ids, random_num = 10,
                             random_length = 100)

randomResults 

width(randomResults)

The function returns a GRangesList object. Every list element has a total length of 100 nt and represents a region on mRNA but possibly separated by intron parts. The name of each element is its relative transcript id.

randomizeTransByOrder function

This is another randomization function. The difference of this function and randomizeTx is that it does not take an argument about transcript ids, instead, it takes a GRangesList region set and picks only one random region within the space indicated by every GRangesList element. So the number of regions in a set this function returns is equal to the number of elements in the input GRangesList object.

Step 1: specify a randomization space

We pick a set of five mRNAs to be the scope of randomization.

trans.ids <- c("170", "782", "974", "1364", "1387")
exons.tx0 <- exonsBy(txdb)
regions.A <- exons.tx0[trans.ids]
regions.A

Step 2: randomizeTransByOrder function

randomizeTransByOrder requires two inputs: a region set with GRangesList format as the randomization scope, and the length of each random region is going to be picked from every GRangesList element.

In the following example, we input a region set containing five regions. The function picks five 100 nt regions from each of them in order. 

set.seed(12345677)
random.regions.A <- randomizeTransByOrder(regions_A = regions.A, 
random_length = 100)
width(regions.A)
width(random.regions.A)

Randomize features into transcriptome

Functions introduced in this section support randomizing two kinds of RNA-related genomic features:

knitr::include_graphics("figures/section4.1.png")

knitr::include_graphics("figures/section4.2.png")

randomizeFeatruesTx function

This function is for randomizing a feature set that does not have isoform ambiguity problem.

Step 1: Specify a feature set (without isoform ambiguity)

We expect users to input features following the format stated in section 3, i.e., either to be a GRanges or GRangesList object with corresponding required metadata.

We generate a region set containing 100 regions of 100nt long as an example feature set.

txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
trans.ids<- c("170", "782", "974", "1364", "1387")
RS1 <-  randomizeTx(txdb, trans.ids, random_num = 100, random_length = 100)

Step 2: randomizeFeaturesTx function

Once we have a feature set, we can use randomizeFeaturesTx function to randomize it. This function takes each feature, randomizes its position only within its corresponding transcript that it belongs to and keeps the original size and other associated data.

randomResults <- randomizeFeaturesTx(RS1, txdb, type = "mature")

By setting the argument type, one can randomize input feature set into different types of transcriptome space. The above example picks random regions from mRNAs. The below example only picks random regions from CDS.

randomResults <- randomizeFeaturesTx(RS1, txdb, type = "CDS")

By setting the argument N, one can get multiple sets of randomized region sets of the input feature set. 

randomResults <- randomizeFeaturesTx(RS1, txdb, N = 5)

randomizeFeatruesTxIA function

This function is for randomizing a feature set that has isoform ambiguity.

Step 1: Specify features with isoform ambiguity

We expect users to input features following the format stated in section 3. For features with isoform ambiguity, users do not need to provide transcript information.

Here we load an m$^6$A dataset as an input feature.

file <- system.file(package="RgnTX", "extdata", "m6A_sites_data.rds")
m6A_sites_data <- readRDS(file)
m6A_sites_data[1:5]

Step 2: randomizeFeaturesTxIA function

One does not know the specific transcript that each m$^6$A feature is associated with. When doing randomization, this function takes each feature and randomizes its position only within the set of isoform transcripts that each feature can be mapped to. We input five features and get five random regions. Each random region is corresponding to each feature.

randomResults <- randomizeFeaturesTxIA(RS = m6A_sites_data[1:5], 
                                       txdb, type = "mature", N = 1)
length(randomResults)

Permutation test function

The most important task of RgnTX is to analyze transcriptome region set association via a permutation test framework. There are several permutation test functions. They require different input ways and support user's different needs. All the functions finally return a report that reflects the level of association being tested and can be visualized by a plot function.

Test between two region sets

This is the most basic strategy. It requires users to input a feature set (to be randomized) and a region set (to be compared with), and evaluates if there is an association between them.

permTestTx function

The permTestTx function is for features without isoform ambiguity. As an example, we generate two region sets A and B. As shown below, by our design they have at least 50 overlaps. We feed them into permTestTx.

set.seed(12345677)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
exons.tx0 <- exonsBy(txdb)
trans.ids <- sample(names(exons.tx0), 50)

A <- randomizeTx(txdb, trans.ids, random_num = 100, random_length = 100)
B <- c(randomizeTx(txdb, trans.ids, random_num = 50, random_length = 100), 
       A[1:50])

permTestTx_results <- permTestTx(RS1 = A, 
                                 RS2 = B, 
                                 txdb = txdb, 
                                 type = "mature",
                                 ntimes = 50, 
                                 ev_function_1 = overlapCountsTx, 
                                 ev_function_2 = overlapCountsTx)

file <- system.file(package="RgnTX", "extdata", "permTestTx_results2.rds")
permTestTx_results <- readRDS(file)

This function involves these arguments:

Among the above arguments, users must at least provide a feature set RS1, a region set RS2 and a TxDb object to the main function.

The permTestTx function automatically randomizes the feature set RS1 and creates a set of randomized region sets, which is outputted as a list named as RSL (by its inner randomization function randomizeFeatruesTx). permTestTx uses the first evaluation function ev_function_1 to calculate the test statistic between RS1 and RS2, and the second evaluation function to compute evaluations between each random region set in RSL and RS2. Finally, all these calculated results are returned as a list, which is defined to be a permTestTx.results object.

names(permTestTx_results)

To get a more intuitive view of the result, we can plot it by plotPermResults function.

p1 <- plotPermResults(permTestTx_results, binwidth = 1)
p1

The number of overlaps between A and B, i.e., orig.ev, should be at least 50. As results show, the random evaluations are much less than orig.ev. A p-value is calculated based on the number that the random evaluations are lower than real observed evaluation, which is less than the significance level 0.05 we set. Additionally, a z-score is computed, which is the distance between mean of random evaluations and real observed evaluation divided by standard deviation of random evaluations. The low p-value and high z-score suggest that there is a statistically significant association between A and B.

permTestTxIA function

Different from permTestTx, this function is for features that have isoform ambiguity, i.e., features that one cannot make sure which specific transcript they are located on. For such feature, usually one can only know a several of isoforms that they can be mapped to. The permTestTxIA function randomizes such features into the corresponding set of isoforms by randomizeFeatruesTxIA function. It conducts the same protocol evaluating an association as function permTestTx does and also returns a permTestTx.results object.

Test involving customPick function

This strategy needs users to provide a feature set and a customPick function. The customPick function should aim to pick a specific kind of regions over a transcript (taking transcript ids as input and outputting a region over each transcript). The aim of this kind of test functions is to evaluate association between input feature set and the customPick regions. For example, suppose one wants to know if some features are associated with the last 200 nt region of transcripts more than expected, one can write a function that picks the transcript’s last 200 nt, and pass this function as a customPick function to the test function, which will then conduct a permutation test to evaluate their association level. Writing this kind of customPick function (specifiy certain regions over a transcript) can be easily done with the help of shiftTx function in section 9.

The input features may be with or without isoform ambiguity. For a feature without isoform ambiguity, default evaluation function (overlapCountsTx) checks if it has overlap with the customPick region (indicated in red) on the transcript where it is located on (see Fig. \@ref(fig:Fig8), a). If the answer is yes, the total number of overlaps is added to one and otherwise zero. While for a feature with isoform ambiguity, default evaluation function (overlapCountsTxIA) maps it with the customPick regions of all its isoforms (see Fig. \@ref(fig:Fig8), b), and every overlap is added not one, but one divided by the number of its isoforms, to the total number of overlaps. In this example, this feature is associated with three isoforms and overlapping with the customPick regions of two of them (indicated in red). Then a value of 2/3 is added to the total overlap.

knitr::include_graphics("figures/section5.2.png")

Function permTestTx_customPick is for features without isoform ambiguity, while function permTestTxIA_customPick is for features with this problem. In this kind of strategies, features are only mapped to the regions in transcripts that they are related to. More details and examples are stated below.

permTestTx_customPick function

As an example, we take a feature set from the CDS space, which contains 100 regions and each region is 200 nt long.

Suppose we do not know how these features are generated. We want to test if these features are associated with the CDS part of transcripts more than expected. To do that, we write a customPick function that can pick the CDS part of these transcripts.

getCDS = function(trans_ids, txdb,...){
  cds.tx0 <- cdsBy(txdb, use.names=FALSE)
  cds.names <- as.character(intersect(names(cds.tx0), trans_ids))
  cds = cds.tx0[cds.names]
  return(cds)
}

The customPick function should involve at least an argument trans_ids and an ellipsis (...) so that other arguments can be passed to the main test function through this ellipsis operator. The permTestTx_customPick function takes the feature set RS1 and the customPick function getCDS.

set.seed(12345677)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
RS1 <- randomizeTx(txdb, "all", random_num = 100, random_length = 200, 
                   type = "CDS")

permTestTx_results <- permTestTx_customPick(RS1 = RS1,
                                            txdb = txdb,
                                            customPick_function = getCDS,
                                            ntimes = 50,
                                            ev_function_1 = overlapCountsTx,
                                            ev_function_2 = overlapCountsTx)
p1 <- plotPermResults(permTestTx_results, binwidth = 1)
p1

file <- system.file(package="RgnTX", "extdata", "permTestTx_results3.rds")
permTestTx_results <- readRDS(file)
p1 <- plotPermResults(permTestTx_results, binwidth = 1)
p1

As Fig. \@ref(fig:Fig9) shows, it is just as we expected that the overlap number between features and the picked CDS is 100, indicated by a red line, which squares with the fact that all these features are picked from the CDS.

permTestTxIA_customPick function

This function receives a feature set (with isoform ambiguity) and a customPick function.

Let’s take an example about testing association between N$^6$-Methyladenosine and the stop codon regions. The m6A sites are derived from an miCLIP-seq dataset (Linder, et al., 2015). Due to technical restriction, it is not clear which specific transcript an m$^6$A feature is located on.

To test this association, we wrote a customPick function getStopCodon picking the stop codon regions of transcripts. It actually picks the last 100 nt of CDS and the first 100 nt of 3'UTR of transcripts given transcript ids. If a given transcript's CDS or 3'UTR part is less than 100 nt, the function will just pick this less than 100 nt region. This function actually uses the function shiftTx in section 9 provided by RgnTX to realize these tasks.

As a quick example, we feed 500 out of 9181 m6A sites and this getStopCodon function into function permTestTxIA_customPick. Its default evaluation function calculates a weighted overlapping number between m6A sites and the stop codon regions.

set.seed(12345677)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
file <- system.file(package="RgnTX", "extdata", "m6A_sites_data.rds")
m6A_sites_data <- readRDS(file)
RS1 <- m6A_sites_data[1:500]

permTestTx_results <- permTestTxIA_customPick(RS1 = RS1,
                                       txdb = txdb,
                                       customPick_function = getStopCodon,
                                       ntimes = 50,
                                       ev_function_1 = overlapCountsTxIA,
                                       ev_function_2 = overlapCountsTx)
p_a <- plotPermResults(permTestTx_results, binwidth = 1)
p_a

txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
file <- system.file(package="RgnTX", "extdata", "permTestTx_results.rds")
permTestTx_results <- readRDS(file)
p_a <- plotPermResults(permTestTx_results, binwidth = 1)
p_a

In addition, if one wants to test the whole 9181 m6A sites with the stop codon, use the following codes, which will of course take a few more minutes than the case involves less data. 

set.seed(12345677)
permTestTx_results <- permTestTxIA_customPick(RS1 = m6A_sites_data,
                                       txdb = txdb,
                                       customPick_function = getStopCodon,
                                       ntimes = 50)
p_b <- plotPermResults(permTestTx_results)
p_b

knitr::include_graphics("figures/section5.5.png")

The weighted number of overlaps between 9181 m6A sites and stop codon regions is 804.93 (calculated by default evaluation function in 6.2), while the mean of random evaluations is 365.28 (calculated by default evaluation function in 6.1). The low p-value shows that the association being tested is statistically significant.

References: Linder, B., et al. (2015) Single-nucleotide-resolution mapping of m6A and m6Am throughout the transcriptome. Nat Methods 12(8): 767-772.

Test with user-provided random sets

permTestTx_customAll function

permTestTx_customAll needs users to provide at least three arguments: two region sets and a list of randomized region sets.

We take an easy example to explain how this function should be used. We want to test the association between two region sets come from the same transcript. We pick two region sets RS1 and RS2 in a transcript that has the id "170", and pick random region sets within a larger space containing 5 transcripts.

set.seed(12345677)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
trans.ids1<- c("170")
RS1 <- randomizeTx(txdb = txdb, trans_ids = trans.ids1,
                   random_num = 20, random_length = 100)
RS2 <- randomizeTx(txdb = txdb, trans_ids = trans.ids1,
                   random_num = 20, random_length = 100)
trans.ids2 <-  c("170", "782", "974", "1364", "1387")
RSL <- randomizeTx(txdb = txdb, trans_ids = trans.ids2,
                   random_num = 20, random_length = 100, N = 50)

We feed these three elements into the permTestTx_customAll function.

permTestTx_results <- permTestTx_customAll(RSL = RSL, RS1 = RS1, RS2 = RS2)
p_3 <- plotPermResults(permTestTx_results)
p_3

The random evaluations indicated by gray boxes are the number of overlaps between RS2 with each region set in RSL. We can see the number of overlaps between RS1 and RS2 is generally larger than these random evaluations (see Fig. \@ref(fig:Fig12)). This accords with the fact that they come from a smaller space so that they are associated with each other more than expected.

Evaluation function

Evaluation functions should define a statistic to be tested between two region sets. RgnTX provides some evaluation functions that can be passed to the main test functions in section 5 via parameters ev_function_1 or ev_function_2. Users can also customize their own evaluation function as examples shown in section 76.3](#custom-evaluation-function). All the evaluation functions require two input region sets following a format in section 4.

overlapCountsTx function

Function overlapCountsTx calculates the number of overlaps between two region sets A and B. It has parameter over_trans deciding whether the overlap is counted over transcriptome or genome. If TRUE, regions in the two input sets that have shared genomic intervals but are located on different transcripts are not considered to have overlap with each other. This function also has a parameter count_once, which decides whether the overlap of multiple regions in the second region set B with a region in A should be counted once or multiple times.

overlapCountsTxIA function

Function overlapCountsTxIA accepts a feature set A (with isoform ambiguity) and a region set B (without isoform ambiguity). Each feature in A related to n_1 isoforms in B and overlapped with n_2 B regions will contribute a value of n_2/ n_1 to the total number of overlaps, which may be a non-integer. This function is used in permTestTxIA_customPick function (see section 6.2.2), calculating a weighted number of overlaps between an m6A feature set (with isoform ambiguity) and a customPick region (without isoform ambiguity).

Custom evaluation function

Users can create their own evaluation function. The custom evaluation function must include parameters of two region sets and an ellipsis (...) so that other parameters can be passed to the main test function. Here are two examples.

Instead of calculating overlapping counts, function overlapWidthTx returns the sum of widths of each overlap, i.e., the total number of overlapping nucleotides between two input region sets.

The distanceTx function calculates the mean of the distance from each region in A to the closest region in B. It contains a user-defined argument, beta, which filters out a corresponding percentage of largest distance values so as to avoid extreme large distances which may make estimating results incomparable.

Parameter setting

In this section, we present some examples to show how different parameter settings influence permutation test results. We take the example in section 6.1.1.

set.seed(12345677)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
exons.tx0 <- exonsBy(txdb)
trans.ids <- sample(names(exons.tx0), 50)

A <- randomizeTx(txdb, trans.ids, random_num = 100, random_length = 100)
B <- c(randomizeTx(txdb, trans.ids, random_num = 50, random_length = 100),
       A[1:50])

permTestTx_results <- permTestTx(RS1 = A,
                                 RS2 = B,
                                 txdb = txdb,
                                 type = "mature",
                                 ntimes = 50,
                                 ev_function_1 = overlapCountsTx,
                                 ev_function_2 = overlapCountsTx)

p_a <- plotPermResults(permTestTx_results, binwidth = 1)

p_b <- plotPermResults(permTestTx_results, test_type = "two-sided")

p_c <- plotPermResults(permTestTx_results, alpha = 0.025)

p_d <- plotPermResults(permTestTx_results, alpha = 0.1)

set.seed(12345677)
permTestTx_results_e <- permTestTx( RS1 = A,
                                    RS2 = B,
                                    txdb = txdb,
                                    type = "mature",
                                    ntimes = 50,
                                    ev_function_1 = overlapWidthTx,
                                    ev_function_2 = overlapWidthTx)
p_e <- plotPermResults(permTestTx_results_e, binwidth = 150)
p_e

set.seed(12345677)
permTestTx_results_f <- permTestTx( RS1 = A,
                                    RS2 = B,
                                    txdb = txdb,
                                    type = "mature",
                                    ntimes = 50,
                                    ev_function_1 = distanceTx,
                                    ev_function_2 = distanceTx, beta = 0.3)
p_f <- plotPermResults(permTestTx_results_f)
p_f

set.seed(12345677)
permTestTx_results_g <- permTestTx(RS1 = A,
                                 RS2 = B,
                                 txdb = txdb,
                                 type = "mature",
                                 ntimes = 500,
                                 ev_function_1 = overlapCountsTx,
                                 ev_function_2 = overlapCountsTx)

p_g <- plotPermResults(permTestTx_results_g, binwidth = 1)
p_g

set.seed(12345677)
permTestTx_results_h <- permTestTx(RS1 = A,
                                 RS2 = B,
                                 txdb = txdb,
                                 type = "mature",
                                 ntimes = 2000,
                                 ev_function_1 = overlapCountsTx,
                                 ev_function_2 = overlapCountsTx)

p_h <- plotPermResults(permTestTx_results_h, binwidth = 1)
p_h

knitr::include_graphics("figures/section7.png")

Shifted Z-scores over transcripts

The value of z-score represents how many standard deviations an observed value is away from the mean of random evaluations. If z-score is close to 0, observed value is about on the mean. If z-score is much larger than 0, observed value is on the positive side and far away from the mean of random evaluations. z-score can be a numerical measurement describing the strength of tested association.

Here we provide a function shiftedZScoreTx which shifts regions of interest from their original positions, applies evaluation function reassessing association and finally calculates new z-scores, what we call shifted z-scores. Calculating shifted z-scores can help us to learn if an association is specially linked to the exact position of features. Function plotShiftedZScoreTx can take the output of shiftedZScoreTx and plot different shifted z-scores when the position of features is changed. If we saw a peak at the original point in a plot of shifted z-scores, we would say that this association is closely related to the specific position of corresponding regions of interest.

To be noticed, different from the localZScore function in r BiocStyle::Biocpkg("regioneR"), shiftZScoreTX in this package provides mode of shifting regions of interest over mRNA space. Users access it by assigning thetype variable a "mature" value.

Let's take the case in section 6.2.2 as an example. In this case, results show that m6A modifications and the stop codon regions tend to have a strong positive association. Let's use shiftedZScoreTx to move stop codons around its original position and calculate corresponding shifted z-scores after each moving. 

set.seed(12345677)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
file <- system.file(package="RgnTX", "extdata", "m6A_sites_data.rds")
m6A_sites_data <- readRDS(file)
RS1 <- m6A_sites_data[1:100]
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
permTestTx_results <- permTestTxIA_customPick(RS1 = RS1,
                                              txdb = txdb,
                                              type = "mature",
                                              customPick_function = getStopCodon,
                                              ntimes = 50)
shiftedZScoreTx_results <- shiftedZScoreTx(permTestTx_results,txdb,
                                           type = 'mature',
                                           window = 500,
                                           step = 50,
                                           ev_function_1 = overlapCountsTxIA)
p1 <- plotShiftedZScoreTx(shiftedZScoreTx_results)
p1

knitr::include_graphics("figures/section9.jpg")

We can see a clear peak in the center of the plot which suggests the tested association is strongly affected by the exact positions of stop codons. The stop codon region tested here is formed by the last 100 nt CDS and the first 100 nt 3'UTR. As shown in Fig. \@ref(fig:Fig14), we can see m6A sites are enriched more on 3'UTR than on CDS, since the peak is deviated to the right (3') direction.

Other useful functions

In this section, we introduce a useful function that can calculate positional shift over transcript regions. It is based on the r BiocStyle::Biocpkg("GenomicRanges ") and r BiocStyle::Biocpkg("IRanges") structure. Most of randomization functions in RgnTX utilize it doing calculation about shifting regions.

shiftTx function

This function accepts a feature set following the format indicated in section 4 and outputs a region set from it. Each output region is from each input feature. Besides this region set, it also receives a number of other parameters:

The direction and strand arguments only receive one value so that this function can only process the same type of direction of shifting and the same strand type of regions at one time.

Example: The following example uses shiftTx to pick the last 200 nt CDS regions in five mRNAs (positive strand). The function starts from the end of CDS (The parameter start takes the largest value of each CDS.) and calculates which coordinates they end after shifting 199 nucleotides to the direction of 5'. This shifting is calculated in mRNA. The impact of introns has been excluded during calculation.

library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
# Five transcripts with positive strand.
trans.ids <- c("170", "782", "974", "1364", "1387")

# Take the CDS part of them.
cds.tx0 <- cdsBy(txdb, use.names = FALSE)
cds.p <- cds.tx0[trans.ids]

# The width of the region from each transcript to be picked is 200.
width <- 200

# The start vector 
start = as.numeric(max(end(cds.p)))

R.cds.last200 <- shiftTx(regions = cds.p, start = start, width = width, 
                         direction = "left", strand = "+")
R.cds.last200

The region set we get is the set of last 200 nt CDS regions of the input five mRNAs. It can be converted to a GRangesList object by the GRanges2GRangesList function. 

R.cds.last200.list <- GRanges2GRangesList(R.cds.last200)
R.cds.last200.list

Each region picked from each feature is 200 nt long.

width(R.cds.last200.list)

This function will not shift or pick regions outside the scope of input features. For example, suppose an input feature is 100 nt. If we assign start a value of the coordinate of its 81th nucleotide and assign width to be 30, the function will only return this feature's last 20 nt region to us.

Performance of randomization function

We evaluate the time it takes for performing randomization functions randomizeTx (Table 1) and randomizeFeaturesTx (Table 2) for one time. We do each group of tests for 10 times and record the average time they spend. We set different size (row) and different number of features (column) for each group. The unit is seconds. The length of randomized regions have slight influence on the randomization time while the number of them has a more direct impact.

table1 <- matrix(c(100, 500, 1000, 2000, 5000, 10000,
         1.783995, 2.027494, 2.341810, 3.081001, 5.097162, 8.434085,
         1.786938, 2.024690, 2.367650, 3.006552, 5.028154, 8.959244,
         1.794495, 2.075795, 2.390498, 3.020262, 5.107760, 9.162599,
         1.761639, 2.046817, 2.585473, 3.302140, 5.187891, 9.148554), nrow = 6, ncol = 5)
table1 = data.frame(table1)
names(table1) <- c('Num/Length', 10,50, 100, 200)
knitr::kable(table1, "simple", caption = 'Randomization time of randomizeTx.')

table2 <- matrix(c(100, 500, 1000, 2000, 5000, 10000,
         1.770159, 2.004670, 2.335991, 3.034632, 5.039495, 8.644615,
         1.801999, 2.024997, 2.322497, 3.042497, 5.194495, 8.500142,
         1.755498, 2.096545, 2.377498, 3.277249, 5.190499, 8.819886,
         1.771499, 2.026999, 2.359494, 3.291889, 5.235691, 9.216726), nrow = 6, ncol = 5)
table2 = data.frame(table2)
names(table2) <- c('Num/Length', 10,50, 100, 200)
knitr::kable(table2, "simple", caption = 'Randomization time of randomizeFeaturesTx.')

Session Info

sessionInfo()


yue-wang-biomath/RgnTX documentation built on Aug. 24, 2023, 1:12 p.m.

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.

Close