YuLab-SMU / ChIPseeker

:dart: ChIP peak Annotation, Comparison and Visualization
https://onlinelibrary.wiley.com/share/author/GYJGUBYCTRMYJFN2JFZZ?target=10.1002/cpz1.585
224 stars 74 forks source link

plot from plotPeakProf2 exon end as central wrongly shows "TSS" not "exon end" #182

Closed HUNNNGRY closed 2 years ago

HUNNNGRY commented 2 years ago

Prerequisites

Describe you issue

Ask in right place

HUNNNGRY commented 2 years ago

here is the code: ChIPseeker::plotPeakProf2(peak = peak, TxDb = txdb, upstream = 1000, downstream = 1000, ignore_strand = F, by = "exon", type = "end_site", nbin = 50,

facet = "row",

                              resample = 500, conf = 0.95

)

MingLi-929 commented 2 years ago

Since I can not get full information from your word, I just give a simple example to explain your question. plotPeakProf2 use tagMatrix() to get the peak count. https://github.com/YuLab-SMU/ChIPseeker/blob/537132283d324bc380d76c654742f0aab73d6056/R/plotTagMatrix.R#L568-L607 and tagMatrix() use getBioRegion() to get the region https://github.com/YuLab-SMU/ChIPseeker/blob/537132283d324bc380d76c654742f0aab73d6056/R/tagMatrix.R#L24-L44

So the output of getBioRegion() can test your question.

files <- ChIPseeker::getSampleFiles()

peak <- files[[4]]
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene::TxDb.Hsapiens.UCSC.hg19.knownGene

region <- ChIPseeker::getBioRegion(TxDb = txdb,
                                   upstream = 1000, 
                                   downstream = 1000,
                                   by = "exon", 
                                   type = "end_site")

GenomicRanges::granges(region)[1]

> GenomicRanges::granges(region)[1]
GRanges object with 1 range and 0 metadata columns:
      seqnames      ranges strand
         <Rle>   <IRanges>  <Rle>
  [1]     chr1 11227-13227      +
  -------
  seqinfo: 93 sequences from an unspecified genome; no seqlengths

Here we get the end region of the first exon. Now we get the exon region to see if it is the end region of the exon

exon <- unlist(GenomicFeatures::exonsBy(txdb))[1]
exon

> exon
GRanges object with 1 range and 3 metadata columns:
    seqnames      ranges strand |   exon_id   exon_name
       <Rle>   <IRanges>  <Rle> | <integer> <character>
  1     chr1 11874-12227      + |         1        <NA>
    exon_rank
    <integer>
  1         1
  -------
  seqinfo: 93 sequences (1 circular) from hg19 genome

So you can see chr1 11227-13227 is the expand region of the end of chr1 11874-12227, which is the region of the first exon.

HUNNNGRY commented 2 years ago

Thank you for your quick reply ! The question is when you further plot with code: ChIPseeker::plotPeakProf2(peak = peak, TxDb = txdb, upstream = 1000, downstream = 1000, ignore_strand = F, by = "exon", type = "end_site", nbin = 50, resample = 500, conf = 0.95 ) seems the central x axis label still shows "TSS" not something like "exon_end", the actual meaning of central point should be the end of exon end not TSS, is that right ?

MingLi-929 commented 2 years ago

It seems that the output of the code is "TTS" rather than "TSS", which means transcription termination site.

files <- ChIPseeker::getSampleFiles()

peak <- files[[4]]
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene::TxDb.Hsapiens.UCSC.hg19.knownGene

> ChIPseeker::plotPeakProf2(peak=peak,
+                           TxDb = txdb,
+                           upstream = 1000, 
+                           downstream = 1000,
+                           ignore_strand = F,
+                           by = "exon", type = "end_site",
+                           nbin = 50,
+                           #facet = "row",
+                           resample = 500, conf = 0.95)
>> binning method is used...2022-03-07 15:47:54
>> preparing end_site regions by exon... 2022-03-07 15:47:54
>> preparing tag matrix by binning...  2022-03-07 15:47:54 
>> Running bootstrapping for tag matrix...       2022-03-07 15:47:56 

exon

HUNNNGRY commented 2 years ago

sorry for the misspelling. So in such situation "TTS" means collapsed end site of all exon ? Also, I'm wondering the plot generated from the code below:

files <- ChIPseeker::getSampleFiles() peak <- files[[4]] txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene::TxDb.Hsapiens.UCSC.hg19.knownGene ChIPseeker::plotPeakProf2(peak = peak, TxDb = txdb, upstream = rel(0.2), downstream = rel(0.2), ignore_strand = F, by = "exon", type = "body", nbin = 50, resample = 500, conf = 0.95 ) "TSS" and "TTS" int this plot mean start and end site of each exon respectively, or they mean transcription start site and transcription end site (and the plot means combine all exons in a transcript without introns) ?

MingLi-929 commented 2 years ago

So in such situation "TTS" means collapsed end site of all exon ?

yes, we get the end site of each exon and expand it with specific number.

"TSS" and "TTS" int this plot mean start and end site of each exon respectively, or they mean transcription start site and transcription end site

They mean start site and end site of each exon. There is stretch of each exon including lengthen or shorten to put each exon with different length in the same plot. So they mean start site and end site of each exon

HUNNNGRY commented 2 years ago

that make sense.
Many thanks for your patient reply