03_plot.R 8.93 KB
Newer Older
1

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

    # Color
44
    n <- length(grlist)
45
46
47
48
49
50
51
52
    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)){
53
        karyoploteR::kpPlotRegions(kp, grlist[[i]], col = colors[i])
54
55
56
    }
    
    # Add legend
57
    if (has_names(grlist)){
58
        graphics::legend('right', fill = colors, legend = names(grlist))
59
60
61
62
63
    }

}


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


107
#' Interval plot GRanges
108
#' @param gr          \code{\link[GenomicRanges]{GRanges-class}}
aditya.bhagwat's avatar
aditya.bhagwat committed
109
110
111
112
#' @param xref        gr var used for scaling x axis
#' @param y           'contig' (default) or name of gr variable
#' @param nperchrom    number (default 1): n head (and n tail) targets 
#'                     shown per chromosome
113
#' @param color_var   'seqnames' (default) or other gr variable
aditya.bhagwat's avatar
aditya.bhagwat committed
114
115
116
#' @param linetype_var NULL (default) or gr variable mapped to linetype
#' @param size_var     NULL (default) or gr variable mapped to size
#' @param facet_var    NULL(default)  or gr variable mapped to facet
aditya.bhagwat's avatar
aditya.bhagwat committed
117
#' @param alpha_var    NULL or gr variable mapped to alpha
aditya.bhagwat's avatar
aditya.bhagwat committed
118
#' @param title        NULL or string: plot title
aditya.bhagwat's avatar
aditya.bhagwat committed
119
#' @param scales       'free', 'fixed', etc
120
121
122
123
124
#' @return ggplot object
#' @seealso  \code{\link{plot_karyogram}}
#' @examples 
#' # SRF sites
#'     require(magrittr)
125
#'     bsgenome <- BSgenome.Mmusculus.UCSC.mm10::BSgenome.Mmusculus.UCSC.mm10
126
#'     bedfile <-  system.file('extdata/SRF.bed',  package = 'multicrispr')
127
128
129
130
131
#'     targets   <- bed_to_granges(bedfile, 'mm10', plot = FALSE)
#'     plot_intervals(targets)
#'     targets %<>% extend(plot = TRUE)
#'     spacers <- find_spacers(targets, bsgenome)
#'     specific <- filter_target_specific(spacers, targets, bsgenome)
aditya.bhagwat's avatar
aditya.bhagwat committed
132
#'     efficient <- filter_efficient(spacers, bsgenome, cutoff = 0.4)
133
134
#'     
#' # PE targets
135
136
137
138
139
140
#'     bsgenome <- BSgenome.Hsapiens.UCSC.hg38::BSgenome.Hsapiens.UCSC.hg38  
#'     gr <- char_to_granges(c(PRNP = 'chr20:4699600:+',
#'                             HBB  = 'chr11:5227002:-',
#'                             HEXA = 'chr15:72346580-72346583:-',
#'                             CFTR = 'chr7:117559593-117559595:+'), 
#'                           bsgenome)
141
#'     spacers <- find_pe_spacers(gr, bsgenome)
142
143
#'     plot_intervals(gr)
#'     plot_intervals(extend_for_pe(gr))
144
#'     plot_intervals(spacers)
145
#'     
146
#' @export
147
plot_intervals <- function(
148
149
150
151
152
153
154
155
156
157
158
    gr, 
    xref         = 'targetname',
    y            = 'names',
    nperchrom    = 1,
    color_var    = 'targetname', 
    facet_var    = 'seqnames', 
    linetype_var = NULL, 
    size_var     = NULL, 
    alpha_var    = NULL,
    title        = NULL, 
    scales       = 'free'
159
){
aditya.bhagwat's avatar
aditya.bhagwat committed
160
161
162
    # Comply
    edge <- targetname <- NULL
    
163
    # Assert, Import, Comply
164
    assert_is_all_of(gr, 'GRanges')
165
    if (!is.null(color_var)) assert_is_a_string(color_var)
166
    assert_is_subset(color_var, names(as.data.table(gr)))
167
168
    contig <- .N <- .SD <- seqnames <- start <- NULL
    strand <- tmp <- width <- xstart <- xend <- . <- NULL
169

170
    # Identify contigs and order on them
171
172
173
    # gr$contig <- GenomicRanges::findOverlaps(
    #                 gr, maxgap = 30, select = 'first', ignore.strand = TRUE)
    # gr %<>% extract(order(.$contig))
174
175
    
    # Prepare plotdt
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
    plotdt <- prepare_plot_intervals(gr, xref, y, nperchrom)
    
    # Core Ranges
    p <-ggplot( plotdt, 
                aes_string(
                    x = 'xstart', xend = 'xend', y = 'y', yend = 'y', 
                    color = color_var, linetype = linetype_var, 
                    size = size_var, alpha = alpha_var)) + 
        facet_wrap(facet_var, scales = scales) + 
        geom_segment(
            arrow = arrow(length = unit(0.1, "inches")))
    
    # Targets
    if (all(c('targetstart', 'targetend') %in% names(mcols(gr)))){
        p <-p + 
            geom_point(aes_string(
                x = 'xtargetstart', y = 'y'), shape = '|', size = 4) +
            geom_point(aes_string(
                x = 'xtargetend',   y = 'y'), shape = '|', size = 4)
    }
    
    # Extensions
    if ('extension' %in% names(mcols(gr))){
        p <-p + geom_segment(aes_string(x='extstart', xend='extend', size = NULL), 
                            linetype = 'dotted',
                arrow = arrow(length = unit(0.1, "inches")))
    }
    p <- p + theme_bw()  +  xlab(NULL)  +  ylab(NULL)  +  ggtitle(title)
    
    # Return
    # print(p)
    p
}

prepare_plot_intervals <- function(gr, xref, y, nperchrom){
    
    # Comply
    edge <- targetname <- xstart <- xend <- width <- NULL
    targetstart <- targetend <- xtargetstart <- xtargetend <- NULL
    extstart <- primer <- revtranscript <- extension <- tmp <- NULL
    
    # Main Ranges
218
    plotdt <- data.table::as.data.table(gr) %>% cbind(names = names(gr))
219
    plotdt %<>% extract(order(seqnames, start))
220
221
222
223
    head_tail <- function(x, n=nperchrom) x [ x %in% c(head(x, n), tail(x, n)) ]
    plotdt %<>% extract(, edge := targetname %in% head_tail(unique(targetname)), by = 'seqnames')
    plotdt %<>% extract(edge==TRUE)
    plotdt %>%  extract(, y      := min(start), by = y)
224
    plotdt %>%  extract(, y      := factor(format(y, big.mark = " ")))
225
    plotdt %>%  extract(, xstart := start-min(start), by = xref)
226
    plotdt %>%  extract(, xend   := xstart + width)
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
    
    # Target marks
    if (all(c('targetstart', 'targetend') %in% names(mcols(gr)))){
        plotdt %>% extract(, xtargetstart := xstart + targetstart-start)
        plotdt %>% extract(, xtargetend   := xend   + targetend-end  )
    }
    
    # Extensions
    if ('extension' %in% names(mcols(gr))){
        plotdt %>% extract(strand=='+', extstart := xstart+18-nchar(primer)[1])
        plotdt %>% extract(strand=='-', extstart := xend-16-(nchar(revtranscript)[1]-1))
        plotdt %>% extract(           , extend   := extstart + nchar(extension)[1]-1)
    }

    # Flip for arrow direction    
242
243
244
    plotdt %>%  extract(strand=='-', tmp    := xend)
    plotdt %>%  extract(strand=='-', xend   := xstart)
    plotdt %>%  extract(strand=='-', xstart := tmp)
245
246
247
248
249
    if ('extension' %in% names(mcols(gr))){
        plotdt %>%  extract(strand=='+', tmp := extend)
        plotdt %>%  extract(strand=='+', extend := extstart)
        plotdt %>%  extract(strand=='+', extstart := tmp)
    }
250
251
    plotdt %>%  extract(, tmp := NULL)
    
252
253
    # Return
    plotdt
254
}