[BioC] Gviz - Visualize Reads - Color by Presence of Junction?

[Hahne, Florian]

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