Gviz - Visualize Reads - Color by Presence of Junction?
1
0
Entering edit mode
Andrew Jaffe ▴ 120
@andrew-jaffe-4820
Last seen 9.6 years ago
Hey, I started playing around with Gviz for visualizing aligned RNAseq reads and have been pretty impressed, particularly with integrating external annotation with read-level data. I've been able to get multi-paneled plots depicting reads and corresponding coverage in specific regions - I wanted to a) confirm that reads from the same fragment are what are connected with lines (and do not indicate junctions, unlike IGV) and b) know if it was possible to color reads that contain junctions in a different color than reads without a junction (e.g. color by the `ngap` column being greater than 0 when `readGAlignmentPairs` or `readGAlignments` is used to read in a bam file) I tried looking online but had a hard time finding anything Thanks, Andrew [[alternative HTML version deleted]]
Coverage Gviz Coverage Gviz • 2.2k views
ADD COMMENT
0
Entering edit mode
Robert Ivanek ▴ 730
@robert-ivanek-5892
Last seen 4 months ago
Switzerland
Hi Andrew, a) If you are using the default .import.bam function from Gviz then the lines connect the reads from the same fragment. However the default function is very simple, it ignores the spliced alignments. If you want to visualize such reads you have to modify the import function something like the one below. This one depends on the function extractAlignmentRangesOnReference from the package GenomicAlignments. b) during the import you can also mark the reads which contains the junctions and later during the plotting you can assign a colour to those reads. Best Robert ##--------------- library(Gviz) require(GenomicAlignments) import.bam2 <- 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]] if (parent.env(environment())[["._trackType"]] == "DataTrack") { res <- if (!as.character(seqnames(selection)[1]) %in% names(sinfo$targets)) { mcols(selection) <- DataFrame(score = 0) selection } else { param <- ScanBamParam(what = c("pos", "qwidth"), which = selection, flag = scanBamFlag(isUnmappedQuery = FALSE)) x <- scanBam(file, param = param)[[1]] cov <- coverage(IRanges(x[["pos"]], width = x[["qwidth"]])) if (length(cov) == 0) { mcols(selection) <- DataFrame(score = 0) selection } else { GRanges(seqnames = seqnames(selection), ranges = IRanges(start = start(cov), end = end(cov)), strand = "*", score = runValue(cov)) } } } else { res <- if (!as.character(seqnames(selection)[1]) %in% names(sinfo$targets)) { mcols(selection) <- DataFrame(id = "NA", group = "NA") selection[0] } else { param <- ScanBamParam(what = c("pos", "qwidth", "strand", "qname", "cigar"), which = selection, flag = scanBamFlag(isUnmappedQuery = FALSE)) x <- scanBam(file, param = param)[[1]] irl <- extractAlignmentRangesOnReference(cigar=x[["cigar"]], pos=x[["pos"]]) irl.length <- sapply(irl, length) GRanges(seqnames = seqnames(selection), ranges = unlist(irl), strand = rep(x[["strand"]], irl.length), id = rep(make.unique(x[["qname"]]), irl.length), group = rep(x[["qname"]], irl.length), feature = rep(ifelse(irl.length>1, "gapped", "other"), irl.length)) } } return(res) } ## test bamFile bamFile <- system.file("extdata/test.bam", package = "Gviz") ## next three lines plot the reads in the default way aTrack4 <- AnnotationTrack(range = bamFile, genome = "hg19", name = "Reads", chromosome = "chr1") plotTracks(list(aTrack4), from = 189990000, to = 190000000) ## aTrack5 <- AnnotationTrack(bamFile, genome = "hg19", chromosome = "chr1", name="Reads", importFunction = import.bam2, stream=T) ## one has to modify the mapping between returned GRanges object ## and the AnnotationTrack ## by setting group = group, lines connect reads from the same fragment ## by setting group = id, lines represents junctions aTrack5 at mapping <- list(id="id", group="group", feature="feature") ## by specifying colour for feature = gapped you can highlight junction ## reads plotTracks(list(aTrack5), from = 189990000, to = 190000000, gapped="red") ##--------------- On 10/12/13 19:48, Andrew Jaffe wrote: > Hey, > > I started playing around with Gviz for visualizing aligned RNAseq reads and > have been pretty impressed, particularly with integrating external > annotation with read-level data. I've been able to get multi- paneled plots > depicting reads and corresponding coverage in specific regions - I wanted > to > > a) confirm that reads from the same fragment are what are connected with > lines (and do not indicate junctions, unlike IGV) and > > b) know if it was possible to color reads that contain junctions in a > different color than reads without a junction (e.g. color by the `ngap` > column being greater than 0 when `readGAlignmentPairs` or `readGAlignments` > is used to read in a bam file) > > I tried looking online but had a hard time finding anything > > Thanks, > Andrew > > [[alternative HTML version deleted]] > > _______________________________________________ > 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 >
ADD COMMENT
0
Entering edit mode
Just to add to this, Robert and I have recently been talking about a somewhat more dedicated AlignedReadsTrack class to produce visualisations similar to the ones people are used to from the IGV browser. This will probably take a couple of days to implement but the request pops up all of the time. We'll try our best to get to this in the new year. (First NY resolution. uh uh?) Florian On 12/11/13 2:39 PM, "Robert Ivanek" <robert.ivanek at="" unibas.ch=""> wrote: >Hi Andrew, > >a) If you are using the default .import.bam function from Gviz then the >lines connect the reads from the same fragment. However the default >function is very simple, it ignores the spliced alignments. If you want >to visualize such reads you have to modify the import function something >like the one below. This one depends on the function >extractAlignmentRangesOnReference from the package GenomicAlignments. > >b) during the import you can also mark the reads which contains the >junctions and later during the plotting you can assign a colour to those >reads. > >Best >Robert > >##--------------- > >library(Gviz) >require(GenomicAlignments) > >import.bam2 <- 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]] > if (parent.env(environment())[["._trackType"]] == "DataTrack") { > res <- if (!as.character(seqnames(selection)[1]) %in% > names(sinfo$targets)) { > mcols(selection) <- DataFrame(score = 0) > selection > } > else { > param <- ScanBamParam(what = c("pos", "qwidth"), > which = selection, > flag = scanBamFlag(isUnmappedQuery = >FALSE)) > x <- scanBam(file, param = param)[[1]] > cov <- coverage(IRanges(x[["pos"]], width = x[["qwidth"]])) > if (length(cov) == 0) { > mcols(selection) <- DataFrame(score = 0) > selection > } > else { > GRanges(seqnames = seqnames(selection), > ranges = IRanges(start = start(cov), > end = end(cov)), strand = "*", score = >runValue(cov)) > } > } > } > else { > res <- if (!as.character(seqnames(selection)[1]) %in% > names(sinfo$targets)) { > mcols(selection) <- DataFrame(id = "NA", group = "NA") > selection[0] > } > else { > param <- ScanBamParam(what = c("pos", "qwidth", "strand", >"qname", "cigar"), > which = selection, flag = >scanBamFlag(isUnmappedQuery = FALSE)) > x <- scanBam(file, param = param)[[1]] > irl <- extractAlignmentRangesOnReference(cigar=x[["cigar"]], >pos=x[["pos"]]) > irl.length <- sapply(irl, length) > GRanges(seqnames = seqnames(selection), > ranges = unlist(irl), > strand = rep(x[["strand"]], irl.length), > id = rep(make.unique(x[["qname"]]), irl.length), > group = rep(x[["qname"]], irl.length), > feature = rep(ifelse(irl.length>1, "gapped", >"other"), irl.length)) > } > } > return(res) >} > >## test bamFile >bamFile <- system.file("extdata/test.bam", package = "Gviz") > >## next three lines plot the reads in the default way >aTrack4 <- AnnotationTrack(range = bamFile, genome = "hg19", > name = "Reads", chromosome = "chr1") >plotTracks(list(aTrack4), from = 189990000, to = 190000000) > >## >aTrack5 <- AnnotationTrack(bamFile, genome = "hg19", > chromosome = "chr1", name="Reads", > importFunction = import.bam2, stream=T) >## one has to modify the mapping between returned GRanges object >## and the AnnotationTrack >## by setting group = group, lines connect reads from the same fragment >## by setting group = id, lines represents junctions >aTrack5 at mapping <- list(id="id", group="group", feature="feature") > >## by specifying colour for feature = gapped you can highlight junction >## reads >plotTracks(list(aTrack5), from = 189990000, to = 190000000, gapped="red") > >##--------------- > > > >On 10/12/13 19:48, Andrew Jaffe wrote: >> Hey, >> >> I started playing around with Gviz for visualizing aligned RNAseq reads >>and >> have been pretty impressed, particularly with integrating external >> annotation with read-level data. I've been able to get multi- paneled >>plots >> depicting reads and corresponding coverage in specific regions - I >>wanted >> to >> >> a) confirm that reads from the same fragment are what are connected with >> lines (and do not indicate junctions, unlike IGV) and >> >> b) know if it was possible to color reads that contain junctions in a >> different color than reads without a junction (e.g. color by the `ngap` >> column being greater than 0 when `readGAlignmentPairs` or >>`readGAlignments` >> is used to read in a bam file) >> >> I tried looking online but had a hard time finding anything >> >> Thanks, >> Andrew >> >> [[alternative HTML version deleted]] >> >> _______________________________________________ >> 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
ADD REPLY

Login before adding your answer.

Traffic: 687 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