add new function and fix bugs #190

MingLi-929 · 2022-05-14T12:46:32Z

update 2022.6.2

(1) combine gr and txdb parameter: users can pass self-made gr to txdb.
(2) combine plotPeakProf3() , plotPeakProf2() and plotPeakProf3() : Users can use plotPeakProf() to call plotPeakProf2() and plotPeakProf3().
The usage lay below.
------------------------------------------------------------------------------->
This pull request do following things:
(1) fix several bugs
(2) perfect the annotation of functions
(3) add test files for getTagMatrix() and plotTagMatrix()
(4) getBioRegion() can support UTR regions(3'UTR + 5'UTR), according to #189
(5) add new function makeBioRegionFromGranges() to support make windows from self-made granges object, according to #189.

getBioRegion() can support UTR regions(3'UTR + 5'UTR)

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

txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
threeUTR <- getBioRegion(TxDb = txdb,
                         by = "3UTR",
                         type = "body")

length(threeUTR)

fiveUTR <- getBioRegion(TxDb = txdb,
                        by = "5UTR",
                        type = "body")

length(fiveUTR)

UTR <- getBioRegion(TxDb = txdb,
                    by = "UTR",
                    type = "body")

> length(UTR) == length(threeUTR) + length(fiveUTR)
[1] TRUE

Now, getBioRegion() support UTR regions

add new function makeBioRegionFromGranges()
According to Feature requests about plotPeakProf2 #189, users will need to investigate insulator or enhancer regions, which can not be attained through txdb object. We add new function makeBioRegionFromGranges() to support users investigate insulator or enhancer regions by using self-made granges object.

These self-made granges object served the same role as txdb object. For example, users can prepare a granges object containing enhancer regions. And makeBioRegionFromGranges() can use these self-made granges object to produce the windows("start_site", "end_site" and "body")

# we consider transcript region as enhancer region
# and make self-made granges object
# they can be the same in the form of granges object
enhancer <- transcripts(txdb)[1:5000,]

## we test three kinds of region, start_site, end_site and body
enhancer_body <- makeBioRegionFromGranges(gr = enhancer,
                                          by = "enhancer",
                                          type = "body")

enhancer_start <- makeBioRegionFromGranges(gr = enhancer,
                                           by = "enhancer",
                                           type = "start_site",
                                           upstream = 1000,
                                           downstream = 1000)

enhancer_end <- makeBioRegionFromGranges(gr = enhancer,
                                         by = "enhancer",
                                         type = "end_site",
                                         upstream = 1000,
                                         downstream = 1000)

Also, we modify the getTagMatrix() to call makeBioRegionFromGranges().

# previous getTagMatrix() can only use txdb object to produce function
# now it support the self-made granges object through gr parameter
peak <- getSampleFiles()[[4]]
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene

# make the window by makeBioRegionFromGranges()
enhancer <- transcripts(txdb)[1:5000,]
  
enhancer_body <- makeBioRegionFromGranges(gr = enhancer,
                                            by = "enhancer",
                                            type = "body")

# input the window made from self-made granges object
mt1 <- getTagMatrix(peak = peak,
                        windows = enhancer_body,
                        weightCol = "V5",
                        nbin = 800)

# without the window parameter
# input the self-made granges object directly
mt2 <- getTagMatrix(peak = peak,
                        gr = enhancer,
                        weightCol = "V5",
                        by = "enhancer",
                        type = "body",
                        upstream = rel(0.2),
                        downstream = rel(0.2),
                        nbin = 800)

# input self-made granges object to make the start_site windows
mt3 <- getTagMatrix(peak = peak,
                      weightCol = "V5",
                      gr = enhancer,
                      by = "gene",
                      type = "start_site",
                      upstream = 1000,
                      downstream = 1000)

plotPeakProf2() is also adjusted to support self-made granges object

plotPeakProf2(peak = peak,
                        by = "gene",
                        type = "body",
                        gr = enhancer)

plotPeakProf2(peak = peak,
                        by = "gene",
                        type = "body",
                        gr = enhancer,
                        upstream = rel(0.2),
                        downstream = rel(0.2),
                        nbin = 800)

yuanlizhanshi · 2022-05-19T02:57:15Z

Sorry to bother you.
Could it be possible to accept two gr object inputs in plotPeakProf2, I want to compare the distribution of H3K4me1 in the enhancer region and non-enhancer region.

MingLi-929 · 2022-05-24T16:08:17Z

According to issue#192, plotPeakProf3() is created to received two or more windows. It can not only received self-made granges object, but also supported to use windows from txdb at the same time. In a word, Users can use self-made granges object to compare with the windows from txdb.

PS: this function do not support to compare region in different types. It should all be
one of "start_site", "end_site" and "body".

the usage lay below:

library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(ChIPseeker)
library(GenomicRanges)

peak <- getSampleFiles()[[4]]
peak_list <- getSampleFiles()[4:5]
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene

## accept two granges object
## compare the peak in enhancer regions and non-enhancer

## self-made enhancer region in the form of granges object
enhancer <- transcripts(txdb)[1:5000,]

## self-made non-enhancer region in the form of granges object
non_enhancer <- unlist(fiveUTRsByTranscript(txdb))[1:5000]
  
gr <- list(enhancer,non_enhancer)

plotPeakProf3(peak = peak,
              conf = 0.95,
              gr = gr,
              by = c("enhancer","non-enhancer"),
              windows_name = c("enhancer","non-enhancer"),
              weightCol = "V5",
              type = "start_site",
              upstream = 1000,
              downstream = 1000)

## users can remove the ribbon plot by omit the conf parameter
plotPeakProf3(peak = peak,
              gr = gr,
              by = c("enhancer","non-enhancer"),
              windows_name = c("enhancer","non-enhancer"),
              weightCol = "V5",
              type = "start_site",
              upstream = 1000,
              downstream = 1000)

## you can also input a list of peaks
plotPeakProf3(peak = peak_list,
              gr = gr,
              conf = 0.95,
              by = c("enhancer","non-enhancer"),
              windows_name = c("enhancer","non-enhancer"),
              weightCol = "V5",
              type = "start_site",
              upstream = 1000,
              downstream = 1000)

## if users want to compare the regions obtained 
## from txbd and self-made granges object
## uses can input gr and txdb in the same time
plotPeakProf3(peak = peak,
              gr = gr,
              conf = 0.95,
              by = c("enhancer","gene"),
              windows_name = c("enhancer","gene"),
              weightCol = "V5",
              type = "start_site",
              upstream = 1000,
              downstream = 1000,
              TxDb = txdb)

## check the body region
plotPeakProf3(peak = peak,
              gr = gr,
              conf = 0.95,
              by = c("enhancer","gene"),
              windows_name = c("enhancer","gene"),
              weightCol = "V5",
              type = "body",
              upstream = 1000,
              downstream = 1000,
              TxDb = txdb,
              nbin = 800)

plotPeakProf3(peak = peak,
              gr = gr,
              conf = 0.95,
              by = c("enhancer","gene"),
              windows_name = c("enhancer","gene"),
              weightCol = "V5",
              type = "body",
              TxDb = txdb,
              nbin = 800)

plotPeakProf3(peak = peak,
              gr = gr,
              conf = 0.95,
              by = c("enhancer","gene"),
              windows_name = c("enhancer","gene"),
              weightCol = "V5",
              type = "body",
              upstream = rel(0.2),
              downstream = rel(0.2),
              TxDb = txdb,
              nbin = 800)

plotPeakProf3(peak = peak_list,
              gr = gr,
              conf = 0.95,
              by = c("enhancer","gene"),
              windows_name = c("enhancer","gene"),
              weightCol = "V5",
              type = "body",
              upstream = rel(0.2),
              downstream = rel(0.2),
              TxDb = txdb,
              nbin = 800)

yuanlizhanshi · 2022-05-25T01:01:33Z

Thank you very much.

Wish you all the best in your scientific research.

GuangchuangYu · 2022-05-26T02:50:38Z

the features are good, but I don't like fun1, fun2, fun3, etc.

We need to design a better API.

For the gr, passing the gr to the txdb parameter is OK, as we did this in annotatePeak.

MingLi-929 · 2022-06-02T13:06:15Z

update usage lay below:

library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(ChIPseeker)
library(GenomicRanges)

peak <- getSampleFiles()[[4]]
peak_list <- getSampleFiles()[4:5]
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene

## accept two granges object
## compare the peak in enhancer regions and non-enhancer

## self-made enhancer region in the form of granges object
enhancer <- transcripts(txdb)[1:5000,]

## self-made non-enhancer region in the form of granges object
non_enhancer <- unlist(fiveUTRsByTranscript(txdb))[1:5000]

plotPeakProf(peak = peak,
             conf = 0.95,
             by = c("enhancer","non-enhancer"),
             windows_name = c("enhancer","non-enhancer"),
             weightCol = "V5",
             type = "start_site",
             upstream = 1000,
             downstream = 1000,
             TxDb = list(enhancer,non_enhancer))


## users can remove the ribbon plot by omit the conf parameter
plotPeakProf(peak = peak,
             TxDb = list(enhancer,non_enhancer),
             by = c("enhancer","non-enhancer"),
             windows_name = c("enhancer","non-enhancer"),
             weightCol = "V5",
             type = "start_site",
             upstream = 1000,
             downstream = 1000)


## you can also input a list of peaks
plotPeakProf(peak = peak_list,
             TxDb = list(enhancer,non_enhancer),
             conf = 0.95,
             by = c("enhancer","non-enhancer"),
             windows_name = c("enhancer","non-enhancer"),
             weightCol = "V5",
             type = "start_site",
             upstream = 1000,
             downstream = 1000)


## if users want to compare the regions obtained 
## from txbd and self-made granges object
## uses can input gr and txdb in the same time
plotPeakProf(peak = peak,
             TxDb = list(enhancer,txdb),
             conf = 0.95,
             by = c("enhancer","gene"),
             windows_name = c("enhancer","gene"),
             weightCol = "V5",
             type = "start_site",
             upstream = 1000,
             downstream = 1000)


## check the body region
plotPeakProf(peak = peak,
             TxDb = list(enhancer,txdb),
             conf = 0.95,
             by = c("enhancer","gene"),
             windows_name = c("enhancer","gene"),
             weightCol = "V5",
             type = "body",
             upstream = 1000,
             downstream = 1000,
             nbin = 800)

plotPeakProf(peak = peak,
             TxDb = list(enhancer,txdb),
             conf = 0.95,
             by = c("enhancer","gene"),
             windows_name = c("enhancer","gene"),
             weightCol = "V5",
             type = "body",
             nbin = 800)


plotPeakProf(peak = peak,
             TxDb = list(enhancer,txdb),
             conf = 0.95,
             by = c("enhancer","gene"),
             windows_name = c("enhancer","gene"),
             weightCol = "V5",
             type = "body",
             upstream = rel(0.2),
             downstream = rel(0.2),
             nbin = 800)


plotPeakProf(peak = peak_list,
             TxDb = list(enhancer,txdb),
             conf = 0.95,
             by = c("enhancer","gene"),
             windows_name = c("enhancer","gene"),
             weightCol = "V5",
             type = "body",
             upstream = rel(0.2),
             downstream = rel(0.2),
             nbin = 800)


plotPeakProf(peak = peak,
             upstream = 1000,
             downstream = 1000,
             by = "gene",
             type = "start_site",
             TxDb = txdb,
             nbin = 800)

# test a list of peak files
plotPeakProf(peak = peak_list,
             upstream = 1000,
             downstream = 1000,
             by = "gene",
             type = "start_site",
             TxDb = txdb,
             nbin = 800)
# test body region
# without extension
plotPeakProf(peak = peak_list,
             by = "gene",
             type = "body",
             TxDb = txdb,
             nbin = 800)

# extend with rel object
plotPeakProf(peak = peak_list,
             by = "gene",
             type = "body",
             TxDb = txdb,
             upstream = rel(0.2),
             downstream = rel(0.2),
             nbin = 800)

# extend with actual number
plotPeakProf(peak = peak_list,
             by = "gene",
             type = "body",
             TxDb = txdb,
             upstream = 1000,
             downstream = 1000,
             nbin = 800)

  

plotPeakProf(peak = peak,
             upstream = 1000,
             downstream = 1000,
             by = "gene",
             type = "start_site",
             TxDb = enhancer,
             nbin = 800)

# test a list of peak files
plotPeakProf(peak = peak_list,
             upstream = 1000,
             downstream = 1000,
             by = "gene",
             type = "start_site",
             TxDb = enhancer,
             nbin = 800)
  
# test body region
# without extension
plotPeakProf(peak = peak_list,
             by = "gene",
             type = "body",
             TxDb = enhancer,
             nbin = 800)
  
# extend with rel object
plotPeakProf(peak = peak_list,
             by = "gene",
             type = "body",
             TxDb = enhancer,
             upstream = rel(0.2),
             downstream = rel(0.2),
             nbin = 800)
  
# extend with actual number
plotPeakProf(peak = peak_list,
             by = "gene",
             type = "body",
             TxDb = enhancer,
             upstream = 1000,
             downstream = 1000,
             nbin = 800)

MingLi-929 added 10 commits May 13, 2022 15:55

optimize the functions

4d8b9c6

supported to get 3UTR and 5UTR at the same time

5ed005a

update annotation of getBioRegion()

91b1fb2

add makeBioRegionFromGranges()

42ad8b4

fix bugs

ad59339

make test file to test getTagMatrix() and related functions

672a65e

fix bugs

beee841

fix bugs in plotTagMatrix

a6135ed

add test files for plotTagMatrix

2d79fd8

update man

cb73086

MingLi-929 mentioned this pull request May 16, 2022

getTagMatrix.binning.internal #191

Closed

MingLi-929 added 4 commits May 23, 2022 23:11

add functions of plotPeakProf3()

04108d2

debug for plotTagMatrix

2f0cfcd

fix bugs

86c5a7e

debugs

dd029b1

MingLi-929 mentioned this pull request May 24, 2022

Feature requests about plotPeakProf2 again... #192

Open

combine function and remove parameters

15b343d

GuangchuangYu merged commit 0e7ba84 into YuLab-SMU:master Jun 3, 2022

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

add new function and fix bugs #190

add new function and fix bugs #190

MingLi-929 commented May 14, 2022 •

edited

Loading

yuanlizhanshi commented May 19, 2022

MingLi-929 commented May 24, 2022

yuanlizhanshi commented May 25, 2022

GuangchuangYu commented May 26, 2022

MingLi-929 commented Jun 2, 2022

add new function and fix bugs #190

add new function and fix bugs #190

Conversation

MingLi-929 commented May 14, 2022 • edited Loading

yuanlizhanshi commented May 19, 2022

MingLi-929 commented May 24, 2022

yuanlizhanshi commented May 25, 2022

GuangchuangYu commented May 26, 2022

MingLi-929 commented Jun 2, 2022

MingLi-929 commented May 14, 2022 •

edited

Loading