03_plot.R 9.25 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
#' @param xref        gr var used for scaling x axis
110
#' @param y           'names' (default) or name of gr variable
aditya.bhagwat's avatar
aditya.bhagwat committed
111
112
#' @param nperchrom    number (default 1): n head (and n tail) targets 
#'                     shown per chromosome
113
#' @param nchrom       number (default 6) of chromosomes shown
114
#' @param color_var   'seqnames' (default) or other gr variable
aditya.bhagwat's avatar
aditya.bhagwat committed
115
116
117
#' @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
118
#' @param alpha_var    NULL or gr variable mapped to alpha
aditya.bhagwat's avatar
aditya.bhagwat committed
119
#' @param title        NULL or string: plot title
aditya.bhagwat's avatar
aditya.bhagwat committed
120
#' @param scales       'free', 'fixed', etc
121
122
123
124
125
#' @return ggplot object
#' @seealso  \code{\link{plot_karyogram}}
#' @examples 
#' # SRF sites
#'     require(magrittr)
126
#'     bsgenome <- BSgenome.Mmusculus.UCSC.mm10::BSgenome.Mmusculus.UCSC.mm10
127
#'     bedfile <-  system.file('extdata/SRF.bed',  package = 'multicrispr')
128
129
130
131
132
#'     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
133
#'     efficient <- filter_efficient(spacers, bsgenome, cutoff = 0.4)
134
135
#'     
#' # PE targets
136
137
138
139
140
141
#'     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)
142
#'     spacers <- find_pe_spacers(gr, bsgenome)
143
144
#'     plot_intervals(gr)
#'     plot_intervals(extend_for_pe(gr))
145
#'     plot_intervals(spacers)
146
#'     
147
#' @export
148
plot_intervals <- function(
149
    gr, xref = 'targetname', y = 'names', nperchrom = 2, nchrom = 4, 
aditya.bhagwat's avatar
aditya.bhagwat committed
150
151
    color_var = 'targetname', facet_var = 'seqnames', linetype_var = NULL, 
    size_var = NULL, alpha_var = NULL, title = NULL, scales= 'free'
152
){
aditya.bhagwat's avatar
aditya.bhagwat committed
153
154
155
    # Comply
    edge <- targetname <- NULL
    
156
    # Assert, Import, Comply
157
    assert_is_all_of(gr, 'GRanges')
158
    if (!is.null(color_var)) assert_is_a_string(color_var)
159
    assert_is_subset(color_var, names(as.data.table(gr)))
160
161
    contig <- .N <- .SD <- seqnames <- start <- NULL
    strand <- tmp <- width <- xstart <- xend <- . <- NULL
162
163

    # Prepare plotdt
164
    plotdt <- prepare_plot_intervals(gr, xref, y, nperchrom, nchrom)
165
166
167
    
    # Core Ranges
    p <-ggplot( plotdt, 
aditya.bhagwat's avatar
aditya.bhagwat committed
168
169
170
                aes_string(x = 'xstart', xend = 'xend', y = 'y', yend = 'y', 
                            color = color_var, linetype = linetype_var, 
                            size = size_var, alpha = alpha_var)) + 
171
        facet_wrap(facet_var, scales = scales) + 
aditya.bhagwat's avatar
aditya.bhagwat committed
172
        geom_segment(arrow = arrow(length = unit(0.1, "inches")))
173
174
175
    
    # Targets
    if (all(c('targetstart', 'targetend') %in% names(mcols(gr)))){
aditya.bhagwat's avatar
aditya.bhagwat committed
176
177
178
179
        p <-p + geom_point(aes_string(
                    x = 'xtargetstart', y = 'y'), shape = '|', size = 4) +
                geom_point(aes_string(
                    x = 'xtargetend',   y = 'y'), shape = '|', size = 4)
180
181
182
183
    }
    
    # Extensions
    if ('extension' %in% names(mcols(gr))){
aditya.bhagwat's avatar
aditya.bhagwat committed
184
185
186
        p <-p + geom_segment(
                    aes_string( x = 'extstart', xend = 'extend', size = NULL), 
                                linetype = 'dotted',
187
188
189
190
191
                arrow = arrow(length = unit(0.1, "inches")))
    }
    p <- p + theme_bw()  +  xlab(NULL)  +  ylab(NULL)  +  ggtitle(title)
    
    # Return
aditya.bhagwat's avatar
aditya.bhagwat committed
192
    p # print(p)
193
194
}

aditya.bhagwat's avatar
aditya.bhagwat committed
195
196
197
198
199
200
201
# Identify contigs and order on them
# gr$contig <- GenomicRanges::findOverlaps(
#                 gr, maxgap = 30, select = 'first', ignore.strand = TRUE)
# gr %<>% extract(order(.$contig))



202
203
204
205
head_tail <- function(x, n){
    idx <- x %in% unique( c(head(x, ceiling(n/2)), tail(x, floor(n/2))))
    x[idx]
}
aditya.bhagwat's avatar
aditya.bhagwat committed
206

207
prepare_plot_intervals <- function(gr, xref, y, nperchrom, nchrom){
208
209
210
211
212
213
    
    # Comply
    edge <- targetname <- xstart <- xend <- width <- NULL
    targetstart <- targetend <- xtargetstart <- xtargetend <- NULL
    extstart <- primer <- revtranscript <- extension <- tmp <- NULL
    
214
    # Prepare data.table. Select chromosomes/targets to plot.
215
    plotdt <- data.table::as.data.table(gr) %>% cbind(names = names(gr))
216
    plotdt %<>% extract(order(seqnames, start))
217
218
219
220
221
222
223
224
225
    plotdt$seqnames %<>% droplevels()
    headtailchroms <- head_tail(levels(plotdt$seqnames), nchrom)
    plotdt %<>% extract(headtailchroms, on = 'seqnames')
    plotdt$seqnames %<>% factor(headtailchroms)
    plotdt %<>% extract( # targets
        , .SD[targetname %in% head_tail(unique(targetname), nperchrom)],
        by = 'seqnames')
    
    # Main ranges
226
    plotdt %>%  extract(, y      := min(start), by = y)
227
    plotdt %>%  extract(, y      := factor(format(y, big.mark = " ")))
228
    plotdt %>%  extract(, xstart := start-min(start), by = xref)
229
    plotdt %>%  extract(, xend   := xstart + width)
230
231
232
233
234
235
236
237
238
239
    
    # 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])
aditya.bhagwat's avatar
aditya.bhagwat committed
240
241
242
        plotdt %>% extract(strand=='-',
                            extstart := xend-16-(nchar(revtranscript)[1]-1))
        plotdt %>% extract(, extend   := extstart + nchar(extension)[1]-1)
243
244
245
    }

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