Question: Gviz color alignments by read strand
0
gravatar for Chris Stubben
13 months ago by
Salt Lake City, Utah
Chris Stubben40 wrote:

Is there a way to use fill.reads to color alignments by read strand?  I would like to plot these on a single track

library(Gviz)
alTrack <- AlignmentsTrack(system.file(package="Gviz", "extdata", 
  "gapped.bam"), isPaired=FALSE)
plotTracks(alTrack, from=2997000, to=2998200, chromosome="chr12", 
  fill.reads=c("red", "blue"))

I tried using the strandedBamImport function from a previous question (Gviz - show plus and minus strand from a bam) but get an Error in .Call2("cigar_ranges" ...  

 

gviz • 387 views
ADD COMMENTlink modified 10 months ago by Robert Ivanek650 • written 13 months ago by Chris Stubben40
Answer: Gviz color alignments by read strand
2
gravatar for Robert Ivanek
10 months ago by
Robert Ivanek650
Switzerland
Robert Ivanek650 wrote:

The way how Gviz handles colouring for AlignmentsTrack is a bit special. It recycles colours specified in the fill.reads along the length of unique read ids (uid column). On top of it, if you plot reads as arrows, it also sorts the reads before plotting according to the strand.

Maybe this piece of code does what you want. It is similar to what you suggested but still retain all plotting functionalities for the AlignmentsTrack class. I tested it on few regions (even with gapped reads) and it seems to work.

library(Gviz)

bf <- system.file(package="Gviz", "extdata", "gapped.bam")

chrom <- "chr12"
from <- 2997500
to <- 2998200

alTrack <- Gviz:::.import.bam.alignments(bf, GRanges(chrom, IRanges(from, to)))
alTrack <- AlignmentsTrack(sort(alTrack))

displayPars(alTrack)$fill.reads <- c(plus="#ECAFAE", minus="#AFAEED")[as.numeric(sort(ranges(alTrack)$readStrand))]

plotTracks(alTrack, from=from, to=to, chromosome=chrom)

Let me know if that works for you.

ADD COMMENTlink written 10 months ago by Robert Ivanek650

That works great, thanks.  I just had to add options(ucscChromosomeNames=FALSE)

ADD REPLYlink modified 10 months ago • written 10 months ago by Chris Stubben40
Answer: Gviz color alignments by read strand
0
gravatar for Chris Stubben
13 months ago by
Salt Lake City, Utah
Chris Stubben40 wrote:

I can color by read strand if I subset a data.frame before creating the alignments track.  

bam1 <- readGAlignments(BamFile(system.file(package="Gviz", "extdata", "gapped.bam")))
bamdf <- data.frame(bam1)
names(bamdf)[1] <- "chromosome"
bamdf <- subset(bamdf, start > 2997000 & end < 2998200)
clrs <- ifelse(bamdf$strand == "+", "#ECAFAE", "#AFAEED")
alTrack <- AlignmentsTrack(bamdf, name="Alignments", fill = clrs)
plotTracks(alTrack)
ADD COMMENTlink written 13 months ago by Chris Stubben40

This only works if ALL reads are plotted in the display.

ADD REPLYlink written 13 months ago by Chris Stubben40

Hi Chris,

what exactly do you want to achieve? Do you want to plot individual reads (coloured by the read strand) or coverage (histogram) for both strands separately?

Robert

ADD REPLYlink written 10 months ago by Robert Ivanek650

I have some Digenome-seq data and want to plot potential cut sites that should have reads stacked on opposite strands (so individual reads colored by read strand).  My code works if there's only a few hundred reads that all fit in the plot.

ADD REPLYlink written 10 months ago by Chris Stubben40
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 448 users visited in the last hour