03_plot.R 6.05 KB
Newer Older
1

aditya.bhagwat's avatar
aditya.bhagwat committed
2
#' Karyo/Interval Plot GRanges(List)
3
#' @param grlist \code{\link[GenomicRanges]{GRanges-class}} or 
4
5
6
#'                    \code{\link[GenomicRanges]{GRangesList-class}}
#' @param title plot title
#' @return list (plot_karyogram) or ggplot (plot_intervals)
7
#' @seealso  \code{\link{plot_intervals}}
8
9
#' @examples 
#' # Plot GRanges
10
11
12
#'     bedfile <-  system.file('extdata/SRF.bed',  package = 'multicrispr')
#'     gr <- bed_to_granges(bedfile, 'mm10', plot = FALSE)
#'     plot_karyogram(gr)
13
14
#'   
#' # Plot GRangesList
15
16
17
#'     flanks  <- left_flank(gr)
#'     grlist <- GenomicRanges::GRangesList(sites = gr, flanks = flanks)
#'     plot_karyogram(grlist)
18
#' @export
aditya.bhagwat's avatar
aditya.bhagwat committed
19
plot_karyogram <- function(
20
    grlist, 
21
    title = unique(genome(grlist))
aditya.bhagwat's avatar
aditya.bhagwat committed
22
){
23
24
    
    # Assert
25
    . <- NULL
26
27
    if (methods::is(grlist, 'GRanges')){
        grlist <- GenomicRanges::GRangesList(grlist)
aditya.bhagwat's avatar
aditya.bhagwat committed
28
    }
29
    assert_is_all_of(grlist, 'GRangesList')
30
    
31
    # Extract relevant chromosomes and order them
32
    chroms <- union(seqlevelsInUse(grlist), standardChromosomes(grlist))
33
    stri_extract <- function(stri, pattern){
34
        stri %>% extract(stri_detect_regex(., pattern)) 
35
36
37
38
39
40
41
    }
    chrs1  <- chroms %>% stri_extract('^(chr)?[0-9]$')    %>% sort()
    chrs2  <- chroms %>% stri_extract('^(chr)?[0-9]{2}$') %>% sort()
    chrsXY <- chroms %>% stri_extract('^(chr)?[XY]$')     %>% sort()
    chrsM  <- chroms %>% stri_extract('^(chr)?MT?$')
    orderedchrs <-  c(chrs1, chrs2, chrsXY, chrsM) %>% 
                    c(sort(setdiff(chroms, .)))
42
    genomeranges <- as(seqinfo(grlist)[orderedchrs], "GRanges")
43
44

    # Color
45
    n <- length(grlist)
46
47
48
49
50
51
52
53
    if (n>0){
        hues <- seq(15, 375, length = n + 1)
        colors  <-  grDevices::hcl(h = hues, l = 65, c = 100)[seq_len(n)]
    }
    
    # Plot
    kp <- karyoploteR::plotKaryotype(genomeranges, main = title)
    for (i in seq_len(n)){
54
        karyoploteR::kpPlotRegions(kp, grlist[[i]], col = colors[i])
55
56
57
    }
    
    # Add legend
58
    if (has_names(grlist)){
59
        graphics::legend('right', fill = colors, legend = names(grlist))
60
61
62
63
64
    }

}


65
plot_tracks <- function(grlist){
66
67
68
    
    group <- . <-  NULL
    
69
    if (methods::is(grlist, 'GRangesList')) gr <- unlist(grlist)
70
71
72
73
    genome  <- unique(genome(seqinfo(gr)))
    assert_is_a_string(genome)
    chrom   <- unique(as.character(seqnames(gr)))[1]
    assert_is_a_string(chrom)
74
75
    
    # Find continuum groups
aditya.bhagwat's avatar
aditya.bhagwat committed
76
77
    gr$group <- GenomicRanges::findOverlaps(
                    gr, maxgap = 1, ignore.strand = TRUE, select = 'first')
78
79
80
81
82
83
    
    # Plot
    coretracks <- list( ideogram = Gviz::IdeogramTrack(
                                        chromosome = chrom, 
                                        genome     = genome), 
                        genomeaxis = Gviz::GenomeAxisTrack())
84
    selectedgr   <- subset(gr, group==1) %>% split(names(.))
85
    annottracks  <- mapply( Gviz::AnnotationTrack, selectedgr, name = names(gr))
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
    Gviz::plotTracks(c(coretracks, annottracks), 
                    background.title = 'gray40', 
                    add = TRUE)

}


to_megabase <- function(y){
    z <- vector('character', length(y))
    
    i <- y>1e6
    z[i] <- paste0(round(y[i]*1e-6), ' M')
    
    i <- y>1e3 & y<=1e6
    z[i] <- paste0(round(y[i]*1e-3), ' K')
    
    i <- y<=1e3
    z[i] <- paste0(round(y[i]), 'b')
aditya.bhagwat's avatar
aditya.bhagwat committed
104
    z %>% magrittr::set_names(names(y))
105
106
107
}


108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
#' Interval plot GRanges
#' @param gr        \code{\link[GenomicRanges]{GRanges-class}}
#' @param color_var string: mcol mapped to color.
#' @param title     plot title
#' @return ggplot object
#' @seealso  \code{\link{plot_karyogram}}
#' @examples 
#' # SRF sites
#'     require(magrittr)
#'     bedfile <-  system.file('extdata/SRF.bed',  package = 'multicrispr')
#'     sites   <- bed_to_granges(bedfile, 'mm10', plot = FALSE)
#'     plot_intervals(sites)
#'     
#'     flanks  <- left_flank(gr, plot = FALSE)
#'     sites$color <- 'sites'
#'     flanks$color <- 'flanks'
#'     plot_intervals(c(sites, flanks))
#'     
#' # PE targets
#'     bs <- BSgenome.Hsapiens.UCSC.hg38::BSgenome.Hsapiens.UCSC.hg38  
128
#'     sites <- GRanges(
129
130
131
132
133
#'                seqnames= c(PRNP = 'chr20:4699600',
#'                            HBB  = 'chr11:5227002',
#'                            HEXA = 'chr15:72346580-72346583',
#'                            CFTR = 'chr7:117559593-117559595'),
#'                strand   = c(PRNP = '+', HBB = '-', HEXA = '-', CFTR = '+'), 
134
#'                seqinfo  = seqinfo(bs))
135
136
#'     sites$color <- names(sites)
#'     plot_intervals(sites)
137
#' @export
138
139
140
141
142
plot_intervals <- function(
    gr, 
    color_var = if ('color' %in% names(mcols(gr))) 'color' else 'seqnames', 
    title = NULL
){
143
    
144
    # Assert, Import, Comply
145
146
    assert_is_all_of(gr, 'GRanges')
    assert_is_subset(color_var, names(as.data.table(gr)))
aditya.bhagwat's avatar
aditya.bhagwat committed
147
148
    contig <- group <- .N <- .SD <- seqnames <- start <- NULL
    strand <- tmp <- width <- xstart <- xend <- y <- NULL
149
150

    # Find adjacent ranges    
aditya.bhagwat's avatar
aditya.bhagwat committed
151
152
    gr$contig <- GenomicRanges::findOverlaps(
                    gr, maxgap = 1, select = 'first', ignore.strand = TRUE)
153
    gr %<>% extract(order(gr$contig))
154
155
    
    # Prepare plotdt
156
    plotdt <- as.data.table(gr)
157
158
159
160
161
162
163
164
165
166
167
168
169
    plotdt <- plotdt[ , .SD[contig %in% c(min(contig), max(contig)) ], 
                        by = c('seqnames')]
    plotdt %<>% extract(order(seqnames, start))
    plotdt %>%  extract(, y      := min(start), by = 'contig')
    plotdt %>%  extract(, y      := factor(round(y*1e-6)))
    plotdt %>%  extract(, xstart := start-min(start), by = 'contig')
    plotdt %>%  extract(, xend   := xstart + width)
    plotdt %>%  extract(strand=='-', tmp    := xend)
    plotdt %>%  extract(strand=='-', xend   := xstart)
    plotdt %>%  extract(strand=='-', xstart := tmp)
    plotdt %>%  extract(, tmp := NULL)
    
    # Plot
170
171
172
173
174
175
176
    p <-ggplot( plotdt, 
                aes_string( x = 'xstart', xend = 'xend', y = 'y', yend = 'y', 
                            color = color_var)) + 
        facet_wrap(~ seqnames, scales = 'free') + 
        geom_segment(arrow = arrow(length = unit(0.1, "inches"))) + 
        theme_bw() + 
        xlab('Bases') + ylab('Megabases') + ggtitle(title)
177
178
179
180
181
182

    # Print and return
    print(p)
    p
    
}