mimbpostproc/12-plot_grohmm.Rmd at master · coregenomics/mimbpostproc · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
# Plot groHMM boundaries

```{r Plot groHMM boundaries, cache = TRUE}
load("hist_chr7.RData")
load("genes.RData")
load("gene.RData")
load("calls.RData")

suppressPackageStartupMessages({
    library(Gviz)
})

## Show 500 kb around the gene of interest.
grange <- granges(gene) + 5e5

## Use DataTrack instead of AlignmentsTrack, because the
## AlignmentsTrack smoothing makes the data hard to see and cannot be
## turned off: https://support.bioconductor.org/p/113298/#113366
##
## Subset histogram to the display range.
hist_ol <- subsetByOverlaps(unlist(hist_chr7), grange)
hist_ol <- split(hist_ol, names(hist_ol))
## Calculate max count for consistent scale across histrogram tracks.
scores <- unlist(sapply(hist_ol, score), use.names = FALSE)
ylim <- c(0, quantile(scores, 0.95))
## Display histogram of reads.
tracks_reads <-
    mapply(function(range, name, grange) {
        DataTrack(range = range,
                  type = "h",
                  name = name,
                  alpha = 0.3,
                  alpha.title = 1L,
                  ylim = ylim)
    }, #
    range = hist_ol,
    name = names(hist_chr7),
    MoreArgs = list(grange = grange))
## Display 50 basepair aggregated HMM de novo gene calls.
tracks_hmms <-
    mapply(function(x, name)
        GeneRegionTrack(x, name = name, ,
                        ## Match Gviz:::.getPlottingFeatures()
                        fill = "#0080ff"),
           x = calls, name = names(calls))
names(tracks_hmms) <- paste(names(calls), "hmm", sep = "-")
## Interleave 2 lists of unequal length
## https://stackoverflow.com/a/56091192
interleave <- function(l1, l2) {
    seq <- seq_len(min(length(l1), length(l2)))
    c(rbind(l1[seq],
            l2[seq]),
      l1[-seq],
      l2[-seq])
}
order <- interleave(names(tracks_reads), names(tracks_hmms))
tracks_reads_hmms <- c(tracks_reads, tracks_hmms)[order]

## Display where the gene is on the chromosome using an ideogram.
chromosome <- as.character(seqnames(gene))
track_ideo <- IdeogramTrack(genome = "hg19", chromosome = chromosome)

## Display annotated genes.
track_genes <- GeneRegionTrack(subset(genes, strand == "+"))

## Plot all the tracks together with axis track scale bars.
plotTracks(c(track_ideo,
             GenomeAxisTrack(),
             track_genes,
             GenomeAxisTrack(scale = 1e5),
             tracks_reads_hmms),
           from = start(grange),
           to = end(grange))
```