Thank you Florian. You just made my day.

Quite note though. When 1st tried using your genomic regions and my bam 
file I got this error:
> plotTracks(dTrack, from = 189990000, to = 190000000, col=c("red", 
> "blue"), groups=c("+", "-"), type="hist", col.histogram=NA)
> Error in validObject(.Object) :
>   invalid class "GRanges" object: 1: NROW(ranges(x)) != length(x)
> invalid class "GRanges" object: 2: NROW(strand(x)) != length(x)
When I repeated with a genomic region for which there are reads in my 
bam it worked like a charm (and fast too!).

Cheers,
António

On 11/15/2013 03:35 PM, Hahne, Florian wrote:
> Hi Antonio,
> the default BAM import function used in the Gviz package is not
> particularly smart. It ignores things like strand, etc. However the
> package was designed with a certain amount of flexibility in mind, and you
> can make a lot of things happening by employing user defined functions. In
> this case, you want to change the BAM import function in order to produce
> a DataTrack with two "samples", the forward and the reverse strand. The
> following code should get you to this, or at least fairly close, and you
> may want to improve on that:
>
> library(Gviz)
> library(Rsamtools)
> strandedBamImport <- function (file, selection) {
>      if (!file.exists(paste(file, "bai", sep = ".")))
>          stop("Unable to find index for BAM file '", file, "'. You can
> build an index using the following command:\n\t",
>               "library(Rsamtools)\n\tindexBam(\"", file, "\")")
>      sinfo <- scanBamHeader(file)[[1]]
>      res <- if (!as.character(seqnames(selection)[1]) %in%
>                 names(sinfo$targets)) {
>          mcols(selection) <- DataFrame(score = 0)
>          selection
>      }else {
>          param <- ScanBamParam(what = c("pos", "qwidth", "strand"),
>                                which = selection, flag =
> scanBamFlag(isUnmappedQuery = FALSE))
>          x <- scanBam(file, param = param)[[1]]
>          gr <- GRanges(strand=x[["strand"]], ranges=IRanges(x[["pos"]],
> width = x[["qwidth"]]), seqnames=seqnames(selection)[1])
>          grs <- split(gr, strand(gr))
>          cov <- lapply(grs[c("+", "-")], function(y) coverage(ranges(y),
> width=end(selection)))
>          pos <- sort(unique(unlist(lapply(cov, function(y) c(start(y),
> end(y))))))
>          if(length(pos)==0){
>              mcols(selection) <- DataFrame(plus=0, minus=0)
>              selection
>          }else{
>              GRanges(seqnames = seqnames(selection)[1],
> ranges=IRanges(start=head(pos, -1), end=tail(pos, -1)),
>                      plus=as.numeric(cov[["+"]][head(pos, -1)]),
>                      minus=-as.numeric(cov[["-"]][head(pos, -1)]))
>          }
>      }
>      return(res)
> }
>
>
> bamFile <- system.file("extdata/test.bam", package = "Gviz")
> dTrack <- DataTrack(range=bamFile, genome="hg19", name="Coverage",
> window=-1, chromosome="chr1", importFunction=strandedBamImport,
> stream=TRUE)
> plotTracks(dTrack, from = 189990000, to = 190000000, col=c("red", "blue"),
> groups=c("+", "-"), type="hist", col.histogram=NA)
>
>
>
> Please note that the coverage for the reverse strand reads is somewhat
> arbitrarily set to negative values in order to mirror the coverage around
> 0. Not perfect, but still usable I hope.
> Florian
>
>
>
>
>
> On 11/14/13 10:02 PM, "António domingues" <amjdomingues@gmail.com> wrote:
>
>> Dear Florian and bioconductor list,
>>
>> I have been using Gviz with quite some success but now I need to
>> generate some plots in which I can represent (by color or direction) the
>> strandness of the read. So far I am generating coverage plots from
>> RNA-seq and ChiP-seq data straight from the bam file using the following
>> methods:
>>
>> ## generates coverage plot
>> bTrackRNAcov <- DataTrack(range = bamRNA, genome = gen, chromosome =
>> chr,  name = "RNA-seq", type = "histogram", col.histogram= "#377EB8",
>> fill="#377EB8")
>>
>> ## plot with individual reads
>>     bTrackRNAcov <- DataTrack(range = bamRNA, genome = gen, chromosome =
>> chr,  name = "RNA-seq")
>>
>> Unfortunately I could not find in either of those an option to
>> differentiate reads in sense and anti-sense strand. In an ideal world I
>> would represent in the 1st plot the coverage for the plus strand above a
>> line and for the minus strand below a line. For the second option
>> coloring the reads in the sense and antisense with different colors.
>> Either option would be fine.
>>
>> I read the manual and searched the confines of the internet, and for the
>> life of me could not find a solution - except separating the sense and
>> anti-sense reads in different files. Is there a more elegant way of
>> doing this?
>>
>> Cheers,
>> António
>>
>>> sessionInfo()
>> R version 3.0.2 (2013-09-25)
>> Platform: x86_64-pc-linux-gnu (64-bit)
>>
>> locale:
>>   [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C LC_TIME=en_US.UTF-8
>> LC_COLLATE=en_US.UTF-8
>>   [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
>> LC_PAPER=en_US.UTF-8       LC_NAME=C
>>   [9] LC_ADDRESS=C               LC_TELEPHONE=C
>> LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
>>
>> attached base packages:
>> [1] parallel  grid      stats     graphics  grDevices utils datasets
>> methods   base
>>
>> other attached packages:
>>   [1] GenomicFeatures_1.14.0 AnnotationDbi_1.24.0 Biobase_2.22.0
>> GenomicRanges_1.14.3   XVector_0.2.0
>>   [6] IRanges_1.20.4         BiocGenerics_0.8.0 Gviz_1.6.0
>> RSvgDevice_0.6.4.3     RColorBrewer_1.0-5
>> [11] data.table_1.8.10      scales_0.2.3 ggplot2_0.9.3.1
>> reshape_0.8.4          plyr_1.8
>>
>> loaded via a namespace (and not attached):
>>   [1] biomaRt_2.18.0      Biostrings_2.30.0   biovizBase_1.10.0
>> bitops_1.0-6        BSgenome_1.30.0     cluster_1.14.4
>>   [7] colorspace_1.2-4    DBI_0.2-7           dichromat_2.0-0
>> digest_0.6.3        gtable_0.1.2        Hmisc_3.12-2
>> [13] labeling_0.2        lattice_0.20-24     latticeExtra_0.6-26
>> MASS_7.3-29         munsell_0.4.2       proto_0.3-10
>> [19] RCurl_1.95-4.1      reshape2_1.2.2      rpart_4.1-3
>> Rsamtools_1.14.1    RSQLite_0.11.4      rtracklayer_1.22.0
>> [25] stats4_3.0.2        stringr_0.6.2       tools_3.0.2
>> XML_3.98-1.1        zlibbioc_1.8.0
>>
>> -- 
>> António Miguel de Jesus Domingues, PhD
>> Postdoctoral researcher
>> Deep Sequencing Group - SFB655
>> Biotechnology Center (Biotec)
>> Technische Universität Dresden
>> Fetscherstraße 105
>> 01307 Dresden
>>
>> Phone: +49 (351) 458 82362
>> Email: antonio.domingues(at)biotec.tu-dresden.de
>> --
>> The Unbearable Lightness of Molecular Biology
>>

-- 
António Miguel de Jesus Domingues, PhD
Postdoctoral researcher
Deep Sequencing Group - SFB655
Biotechnology Center (Biotec)
Technische Universität Dresden
Fetscherstraße 105
01307 Dresden

Phone: +49 (351) 458 82362
Email: antonio.domingues(at)biotec.tu-dresden.de
--
The Unbearable Lightness of Molecular Biology


	[[alternative HTML version deleted]]

