Gviz color alignments by read strand
2
0
Entering edit mode
@chris-stubben-12524
Last seen 18 months ago
Salt Lake City, Utah

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 • 911 views
ADD COMMENT
2
Entering edit mode
Robert Ivanek ▴ 700
@robert-ivanek-5892
Last seen 16 hours ago
Switzerland

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 COMMENT
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode
@chris-stubben-12524
Last seen 18 months ago
Salt Lake City, Utah

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 COMMENT
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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 REPLY
0
Entering edit mode

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 REPLY

Login before adding your answer.

Traffic: 373 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6