[BioC] Gviz - Plot genes and data from - strand (3'-5') in 5'-3' direction

Thu Sep 12 17:42:50 CEST 2013

Hi Dominik,

Actually this function assumes that reads object is of DataTrack class
and genes object of AnnotationTrack class.

I did not tested the function a lot, this was just suggestion how to
achieve the plot you were asking for. My understating was that you would
like to flip only count/reads/data which overlap with genes on the minus
strand. Instead of replacing start values in the DataTrack it is
necessary to replace the the entire range. Than you should not get
negative width values:

reverseCounts <- function(reads, genes) {
    genes <- genes[strand(genes)=="-"]
    strand(genes) <- "*"
    ##
    ov <- findOverlaps(ranges(reads), ranges(genes))
    ww <- width(reads)[queryHits(ov)]
    ranges(ranges(reads))[queryHits(ov)] <-
       IRanges(start(genes)[subjectHits(ov)] +
               (end(genes)[subjectHits(ov)] -
                 start(reads)[queryHits(ov)]), width=ww)
    ##
    return(reads)
}

Best
Robert

On 12/09/13 17:35, Jager, Dominik wrote:
> Thanks a lot,
> 
> I'm an absolute bioconductor beginner and I neither have a "bioinformatics" background. So sorry if I'm a little slow ;)
> 
> So can I simply use my GRanges  objects for the annotation and my data with your function? I' m also not really sure, if I understand the function... Also with your "new" function I get this error: dT2 <- reverseCounts(aT, dT1)
> 
> Fehler in validObject(x) : 
>   invalid class “IRanges” object: 'widths(x)' cannot contain negative values
> 
> Dominik
> 
> Dr. Dominik Jäger
> Biological Sciences Department of Microbiology
> 405 Biological Science Building | 484 West 12th Avenue Columbus, OH 43210
> (614)-682-6890 Office
> jager.9 at osu.edu osu.edu
> 
> ________________________________________
> Von: bioconductor-bounces at r-project.org [bioconductor-bounces at r-project.org]" im Auftrag von "Robert Ivanek [robert.ivanek at unibas.ch]
> Gesendet: Donnerstag, 12. September 2013 11:21
> An: bioconductor at r-project.org
> Betreff: Re: [BioC] Gviz - Plot genes and data from - strand (3'-5') in 5'-3' direction
> 
> I realized it is better to change the coordinates than the actual values:
> 
> reverseCounts <- function(reads, genes) {
>     genes <- genes[strand(genes)=="-"]
>     strand(genes) <- "*"
>     ##
>     ov <- findOverlaps(ranges(reads), ranges(genes))
>     ww <- width(reads)[queryHits(ov)]
>     start(reads)[queryHits(ov)] <- start(genes)[subjectHits(ov)] +
>        (end(genes)[subjectHits(ov)] - start(reads)[queryHits(ov)])
>     width(reads)[queryHits(ov)] <- ww
>     ##
>     return(reads)
> }
> 
> Best
> Robert
> 
> 
> On 12/09/13 17:01, Robert Ivanek wrote:
>> I see, sorry I got it wrong.
>> It cannot be done automatically in Gviz but maybe following can help you:
>>
>> library(GenomicRanges)
>> library(Gviz)
>>
>> genes.gr <- GRanges("chr1", IRanges(c(10,20,30,40), width=8),
>>                     strand=c("+","-","+","-"))
>> aT <- AnnotationTrack(genes.gr)
>>
>> reads.gr <- GRanges("chr1", IRanges(seq(1,50,2), width=1),
>>                     strand="+", count=seq(1,50,2))
>> dT1 <- DataTrack(reads.gr)
>>
>> ## this function expects DataTrack in which you would like to flip the
>> points and AnnotationTrack with genes which would define regions you
>> would like to flip, those regions (genes) have to be on minus strand
>>
>> reverseCounts <- function(reads, genes) {
>>     genes <- genes[strand(genes)=="-"]
>>     strand(genes) <- "*"
>>     ##
>>     ov <- findOverlaps(ranges(reads), ranges(genes))
>>     indx <- split(queryHits(ov), subjectHits(ov))
>>     indx <- lapply(indx, rev)
>>     indx <- unlist(indx)
>>     values(reads)[,queryHits(ov)] <- values(reads)[,indx]
>>     ##
>>     return(reads)
>> }
>> ##
>> dT2 <- reverseCounts(dT1, aT)
>> ##
>> plotTracks(list(aT, dT1, dT2))
>>
>> Is that what you want?
>> Best
>> Robert
>>
>>
>> On 12/09/13 16:21, Jager, Dominik wrote:
>>> ...yes, thats correct.
>>>
>>> So what I really want is to plot my data (see attached PNG) flipped by 180 degrees. So now the genes and the data were plotted in 3' to 5' direction, but I want in 5' to 3'/
>>>
>>> Here some code:
>>>
>>>> setwd("~/WORKING/")
>>>> library(Gviz)
>>>> options(ucscChromosomeNames = FALSE)
>>>> tkgff <- import.gff3("~/WORKING/NC_006624.gff", format = "gff3", genome= "NC_06624.1", asRangedData = FALSE)
>>>> atrack <- AnnotationTrack(tkgff, strand = "*", chromosome = chr, genome = gen, name = TK)
>>>> gtrack <- GenomeAxisTrack()
>>>> WIG <- import.wig("~/S0_20min_reverse.wig", asRangedData = FALSE)
>>>> dtrack <- DataTrack(WIG, name = "S0+")
>>> >from <- 1049512
>>>> to <- 1052176
>>>> plotTracks(list(gtrack,atrack, dtrack), from = from, to = to, cex = 1, type = "h")
>>>
>>> I also attached the plot as PNG.
>>>
>>>
>>> Dominik
>>>
>>> Dr. Dominik Jäger
>>> Biological Sciences Department of Microbiology
>>> 405 Biological Science Building | 484 West 12th Avenue Columbus, OH 43210
>>> (614)-682-6890 Office
>>> jager.9 at osu.edu osu.edu
>>>
>>> ________________________________________
>>> Von: bioconductor-bounces at r-project.org [bioconductor-bounces at r-project.org]" im Auftrag von "Steve Lianoglou [lianoglou.steve at gene.com]
>>> Gesendet: Donnerstag, 12. September 2013 09:57
>>> An: Robert Ivanek
>>> Cc: bioconductor at r-project.org list
>>> Betreff: Re: [BioC] Gviz - Plot genes and data from - strand (3'-5') in 5'-3' direction
>>>
>>> Hi,
>>>
>>> On Thu, Sep 12, 2013 at 6:14 AM, Robert Ivanek <robert.ivanek at unibas.ch> wrote:
>>>> Hi Dominik,
>>>>
>>>> You can achieve that by changing the strand in corresponding track
>>>> objects.
>>> [snip]
>>>
>>> While this might change the strand of a gene, I do not think this is
>>> what the OP is really after -- as I have previously been keen to try
>>> to implement the approach that (I think) the OP wants.
>>>
>>> I believe the OP would like to plot "the mirror image" of the plot
>>> that Gviz would "normally" show given the start/stop coordinates that
>>> the plotting functionality infers (or is given).
>>>
>>> The approach you suggest wouldn't properly "flip" my gene models and
>>> the data associated with them -- it would only show them on the
>>> positive strand. The 5'utr of the gene that was originally on the
>>> negative strand would look as if it were the 3'utr of the gene on the
>>> (now annotated) positive strand, and vice versa, and if I were (say)
>>> plotting ChIP-seq data over this locus, the TF would appear to be
>>> binding the 3'utr of a gene on the positive strand (or downstream of
>>> it), when that's not really the case.
>>>
>>> -steve
>>>
>>> --
>>> Steve Lianoglou
>>> Computational Biologist
>>> Bioinformatics and Computational Biology
>>> Genentech
>>>
>>> _______________________________________________
>>> Bioconductor mailing list
>>> Bioconductor at r-project.org
>>> https://stat.ethz.ch/mailman/listinfo/bioconductor
>>> Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor
>>>
>>>
>>> _______________________________________________
>>> Bioconductor mailing list
>>> Bioconductor at r-project.org
>>> https://stat.ethz.ch/mailman/listinfo/bioconductor
>>> Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor
>>>
>>
>> _______________________________________________
>> Bioconductor mailing list
>> Bioconductor at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/bioconductor
>> Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor
>>
> 
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at r-project.org
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor
> 
>