Entering edit mode
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 <- AnnotationTrackgenes.gr)
>>
>> reads.gr <- GRanges("chr1", IRanges(seq(1,50,2), width=1),
>> strand="+", count=seq(1,50,2))
>> dT1 <- DataTrackreads.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
>
>