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
>
