R/plot_genome.R
In shuyuzheng/Chloroplot: A package to plot Chloroplast genome.
Defines functions
PlotMitGenome
PlotPlastidGenome
PlotTab
Documented in
PlotMitGenome
PlotPlastidGenome
PlotTab
#' Generate table for PlotGenome function
#'
#' @param gbfile A character. It is the GI or GenBank accession for the
#' interested specie's chloroplast genome files.
#' @param local.file A logical value. If it is \code{TRUE}, the value passed to
#' \code{gbfile} should be a local GB file.
#' @param gc.window A integer. It indicates the window sized for plot
#' GC count.
#'
#' @return A list. It contains 3 tables:
#' \itemize{
#' \item ir_table: a data frame contains information for IR region.
#' \item gene_table: a data frame contains information of genes.
#' \item gc_count: a data frame contains GC count information.
#' }
#'
#' @export
#' @importFrom magrittr %>%
#' @import dplyr
PlotTab <- function(gbfile, local.file = FALSE, gc.window = 100){
my_env <- new.env(parent = emptyenv())
my_env$gb <- NULL
my_env$def <- NULL
my_env$sp_name <- NULL
my_env$genome <- NULL
my_env$l <- NULL
my_env$gene_table <- NULL
my_env$ir <- NULL
if (local.file){
tryCatch({
my_env$gb <- genbankr::readGenBank(gbfile)
my_env$genome <- genbankr::getSeq(my_env$gb)[[1]]
#genome <- rdnFixer(genome)
my_env$l<- Biostrings::nchar(my_env$genome)
my_env$def <- my_env$gb@definition
my_env$sp_name <- sp.name(my_env$def)
my_env$gene_table <- geneTableRead(my_env$gb, my_env$genome)
}, error = function(e){
my_env$gb <- genbankr::parseGenBank(gbfile)
my_env$def <- my_env$gb$DEFINITION
my_env$sp_name <- sp.name(my_env$def)
my_env$genome <- my_env$gb$ORIGIN[[1]]
my_env$l <- Biostrings::nchar(my_env$genome)
my_env$gene_table <- geneTableParsed(my_env$gb, my_env$genome)
})
} else {
tryCatch({
my_env$gb <- genbankr::readGenBank(text = fetch.gb(gbfile))
my_env$genome <- genbankr::getSeq(my_env$gb)
my_env$l <- Biostrings::nchar(my_env$genome)
my_env$def <- my_env$gb@definition
my_env$sp_name <- sp.name(my_env$df)
my_env$gene_table <- geneTableRead(my_env$gb, my_env$genome)
}, error = function(e){
my_env$gb <- genbankr::parseGenBank(text = fetch.gb(gbfile))
my_env$def <- my_env$gb$DEFINITION
my_env$sp_name <- sp.name(my_env$def)
my_env$genome <- my_env$gb$ORIGIN[[1]]
my_env$l <- Biostrings::nchar(my_env$genome)
my_env$gene_table <- geneTableParsed(my_env$gb, my_env$genome)
})
}
# 1. Organelle ------------------------------------------------------------
plastid <- TRUE
if (grepl("mitochondrion", my_env$def)){
plastid <- FALSE
}
# 2. IR LSC SSC -----------------------------------------------------------
# when seed.size is too short it will raise some error in irDetect function
# If it is too long, it cannot correctly detect the IR resions has mutation
# closing the ends.
if (plastid){
tryCatch({
my_env$ir <- irDetect(my_env$genome, seed.size = 100)
}, error = function(e){
my_env$ir <- irDetect(my_env$genome, seed.size = 1000)
}, warning = function(w) {
my_env$ir <- irDetect(my_env$genome, seed.size = 1000)
})
gene_class <- c("psa","psb","pet","atp","ndh","rbc","rpo","rps","rpl",
"clp|mat|inf","ycf","trn","rrn")
} else {
my_env$ir <- list(ir_table = NULL, indel_table = NULL)
gene_class <- c("nad|nd","sdh","cob","cox|cytb","atp","ccmF","mtt","rps","rpl",
"mat","orf","trn","rrn")
}
my_env$gene_table$class <- rep("OTHER", nrow(my_env$gene_table))
for (i in gene_class){
my_env$gene_table$class[which(grepl(as.character(i),
my_env$gene_table$gene, perl = TRUE,
ignore.case = TRUE))] <- i
}
# 3. GC count -------------------------------------------------------------
if (my_env$l > 500000){
gc.window <- 200
} else if (my_env$l < 100000){
gc.window <- 50
}
gc_count_list <- gc_count(my_env$genome, view.width = gc.window)
gc_count <- gc_count_list[[1]]
gc_total <- gc_count_list[[2]]
gc_count$chr <- rep("chr1", nrow(gc_count))
gc_count <- select(gc_count, chr, position, gc_count)
tables <- list(ir_table = my_env$ir$ir_table, indel_table = my_env$ir$indel_table,
gc_count = gc_count, gc_total = gc_total, plastid = plastid,
sp_name = my_env$sp_name, genome_len = my_env$l, gene_table = my_env$gene_table,
genome = my_env$genome)
return(tables)
}
#' Generate Plasit Genome Plot
#'
#' @param plot.tables a list contains information of IR region, gene, and gc
#' count of the genome. It can be generated by function \code{\link{PlotTab}}.
#' @param save A logical value. If it is \code{TRUE}, the plot will be saved in
#' work directory.
#' @param file.type A charactor. It indicates the format of the file in which
#' the plot will be saved. Options are: "pdf", "png", "jpeg","bmp", "tiff".
#' By default, it is set as "pdf".
#' @param file.name A charactor. It indicates the name of the file in which
#' the plot will be saved. By default, it is set as the specie's name.
#' @param shadow A logical value. If it is \code{TRUE}, the section of IR
#' regions will be highlighted by shadows.
#' @param shadow.color An R color object. It indicates the color for the shadow
#' casted from IR sectors.
#' @param legend A logical value. If it is \code{TRUE}, the legend for gene
#' colors will be shown.
#' @param show.indel A logical value. If it is \code{TRUE}, the SNP, insertion
#' and deletion areas in IR regions will be highlighted on the plot.
#' @param ssc.converse A logical value. If it is \code{TRUE}, the SSC region
#' will be converted to its reverse complementary version
#' @param lsc.converse A logical value. If it is \code{TRUE}, the LSC region
#' will be converted to its reverse complementary version
#' @param ira.converse A logical value. If it is \code{TRUE}, the IRA region
#' will be converted to its reverse complementary version
#' @param irb.converse A logical value. If it is \code{TRUE}, the IRB region
#' will be converted to its reverse complementary version
#' @param background An R color object. It indicates the color for the background of
#' entire plot area.
#' @param height A vector of numeric value. The elements of it indicat the
#' height of gene layer, GC count layer and IR region layer, respectively.
#' Default setting is "0.1, 0.2, 0.07". The taltal circle plot region always
#' has a radius of 1, so a height of 0.1 means 10\% of the circle radius.
#' @param gc.color An R color object. It indicates the color for the lines in gc count
#' plot.
#' @param gc.background An R color object. It indicates the color for the background
#' of gc count plot area.
#' @param info.background An R color object. It indicates the color for the background
#' of central area where species' information was shown.
#' @param ir.color An R color object. It indicates the color for the background
#' of IR sectors.
#' @param ssc.color An R color object. It indicates the color for the background
#' of SSC sectors.
#' @param lsc.color An R color object. It indicates the color for the background
#' of LSC sectors.
#' @param ir.gc A logical value. If it is \code{TRUE}, the gc content for
#' IR, LSC, and SSC regions will be plot with deep colors.
#' @param gc.per.gene A logical value. If it is \code{TRUE}, the gc content for
#' each gene will be plot on the gene layer with deep colors.
#' @param pseudo A logical value. If it is \code{TRUE}, the pseudo genes (if
#' there is some in the species) will be marked with a "*" at the end of the
#' labels.
#' @param cu.bias A logical value. If it is \code{TRUE}, the condon usage bias
#' metric for each gene will be shown in the labels. The metric “Measure
#' Independent of Length and Composition (MILC)” was used for evaluate the bias.
#' You can find more details in the reference paper.
#' @param text.size A numeric value. It indicates the size of all texts in the
#' plot. For exmple, \code{text.size = 1.5} means all the texts in the plot will
#' be enlarged as 1.5 times of their original size.
#' @param gc.per.gene A logical value. If it is \code{TRUE}, the GC content of
#' each gene will be show by a darker part in gene rectangles.
#' @param genome.length A logical value. If it is \code{TRUE}, the length of
#' genome will be shown in the center of the plot.
#' @param total.gc A logical value. If it is \code{TRUE}, the GC content of
#' whole genome will be shown in the center of the plot.
#' @param gene.no A logical value. If it is \code{TRUE}, the number of genes in
#' the genome will be shown in the center of the plot.
#' @param rrn.no A logical value. If it is \code{TRUE}, the number of rRNAs in
#' the genome will be shown in the center of the plot.
#' @param trn.no A logical value. If it is \code{TRUE}, the number of tRNAs in
#' the genome will be shown in the center of the plot.
#' @param organelle A logical value. If it is \code{TRUE}, the organelle type of
#' the genome will be shown in the center of the plot.
#' @param psa.color An R color object. It indicates the color for genes of
#' photosystem I
#' @param psb.color An R color object. It indicates the color for genes of
#' photosystem II
#' @param pet.color An R color object. It indicates the color for genes of
#' cytochrome b/f complex
#' @param atp.color An R color object. It indicates the color for genes of
#' ATP synthesis
#' @param ndh.color An R color object. It indicates the color for genes of
#' NADH dehydrogenase
#' @param rbc.color An R color object. It indicates the color for genes of
#' RubisCO larg subunit
#' @param rpo.color An R color object. It indicates the color for genes of
#' RNA polymerase
#' @param rsp.color An R color object. It indicates the color for genes of
#' small ribosomal protein
#' @param rpl.color An R color object. It indicates the color for genes of
#' large ribosomal protein
#' @param clp_mat_inf.color An R color object. It indicates the color for genes
#' of clpP, matK, or infA
#' @param ycf.color An R color object. It indicates the color for genes of
#' hypothetical reading frame
#' @param trn.color An R color object. It indicates the color for genes of
#' transfer RNA
#' @param rrn.color An R color object. It indicates the color for genes of
#' ribosomal RNA
#' @param other_gene.color An R color object. It indicates the color for other
#' genes
#' @param show.gene A vector of characters. It indicates which classes of genes
#' will be shown on the plot. A valiable values are "psa","psb","pet","atp",
#' "ndh","rbc", "rpo","rps","rpl", "clp|mat|inf","ycf", "trn","rrn", "OTHER"
#' @param gene_axis_ir.color An R color object. It indicates the color for other
#' genes.
#' @param gene_axis_lsc.color An R color object. It indicates the color for other
#' genes.
#' @param gene_axis_ssc.color An R color object. It indicates the color for other
#' genes.
#' @param customize.ring1 A data frame. It must contain 2 columns:
#' \itemize{
#' \item \strong{position}: 1-base genomic coordinate for the features.
#' \item \strong{value}: the values for the features.
#' }
#'
#' @param customize.ring1.type A charactor. It indicate the plot type in
#' customize.ring1. Avaliable values are "line", "line + filling", "line + dot",
#' "line + dot + filling", "step line", "step line + filling", "vertical line"
#' @param customize.ring1.color An R color object. It indicates the color for
#' the plots in customized ring 1.
#' @param customize.ring2 A data frame. It must contain 2 columns:
#' \itemize{
#' \item \strong{position}: 1-base genomic coordinate for the features.
#' \item \strong{value}: the values for the features.
#' }
#' @param customize.ring2.type A charactor. It indicate the plot type in
#' customize.ring2. Avaliable values are "line", "line + filling", "line + dot",
#' "line + dot + filling", "step line", "step line + filling", "vertical line"
#' @param customize.ring2.color An R color object. It indicates the color for
#' the plots in customized ring 2.
#' @param customize.ring3 A data frame. It must contain 2 columns:
#' \itemize{
#' \item \strong{start}: 1-base genomic coordinate for the start point of the
#' features.
#' \item \strong{end}: 1-base genomic coordinate for the end point of the
#' features.
#' \item \strong{value}: the values for the features.
#' }
#' @param customize.ring3.color An R color object. It indicates the color for
#' the plots in customized ring 3.
#' @return A plot for chloroplast genome.
#'
#' @references Supek, Fran, and Kristian Vlahovicek. “Comparison of codon usage
#' measures and their applicability in prediction of microbial gene
#' expressivity.” BMC bioinformatics vol. 6 182. 19 Jul.
#' 2005, doi:10.1186/1471-2105-6-182
#'
#' @importFrom magrittr %>%
#' @import dplyr
#' @export
PlotPlastidGenome <- function(plot.tables, save = TRUE, file.type = "pdf",
text.size = 1, height = c(0.1, 0.2, 0.07),
file.name = NULL, shadow = TRUE, ir.gc = TRUE,
gc.per.gene = TRUE, pseudo = TRUE, legend = TRUE,
ssc.converse = FALSE, lsc.converse = FALSE,
ira.converse = FALSE, irb.converse = FALSE,
genome.length = TRUE, total.gc = TRUE, show.indel = TRUE,
gene.no = TRUE, rrn.no = TRUE,trn.no = TRUE,
organelle_type = TRUE,
background = "grey90",gc.color = "grey30",
gc.background = "grey70", info.background = "black",
ir.color = "#2F3941", ssc.color = "#82B6E2",
lsc.color = "#299E96", shadow.color = "#0000FF20",
psa.color = "#2A6332", psb.color = "#4C8805",
pet.color = "#7F992C", atp.color = "#9FBB3D",
ndh.color = "#FEEE50", rbc.color = "#4D9E3F",
rpo.color = "#AE2D29", rsp.color = "#D6AD7C",
rpl.color = "#9C7A4B", clp_mat_inf.color = "#D9662D",
ycf.color = "#71B8A9", trn.color = "#172C7F",
rrn.color = "#D1382A", other_gene.color = "#7D7D7D",
show.genes = c("psa","psb","pet","atp","ndh","rbc",
"rpo","rps","rpl", "clp|mat|inf","ycf",
"trn","rrn", "OTHER"),
gene_axis_ir.color = ir.color,
gene_axis_ssc.color = ssc.color,
gene_axis_lsc.color = lsc.color,
cu.bias = TRUE, customize.ring1 = NULL,
customize.ring1.type = "line",
customize.ring1.color = gc.color,
customize.ring2 = NULL,
customize.ring2.type = "line",
customize.ring2.color = gc.color,
customize.ring3 = NULL,
customize.ring3.color = gc.color){
# Unpack plot.table
genome <- plot.tables$genome
ir_table <- plot.tables$ir_table
gc_count <- plot.tables$gc_count
gc_total <- plot.tables$gc_total
sp_name <- plot.tables$sp_name
gene_table <- plot.tables$gene_table
genome <- plot.tables$genome
indel_table <- plot.tables$indel_table
l <- Biostrings::nchar(genome)
# 1. Modify gene table
# ssc covert
if (ssc.converse){
if (nrow(ir_table) == 1){
warning("Didn't get IR region from thid species. ",
"It's impossible to convert SSC region.")
} else{
tmp <- convert_region(ir_table = ir_table, l = l, region = "SSC",
genome = genome, gene_table = gene_table,
customize.ring1 = customize.ring1,
customize.ring2 = customize.ring2,
customize.ring3 = customize.ring3)
gene_table <- tmp$gene_table
gc_count <- tmp$gc_count
gc_total <- tmp$gc_total
customize.ring1 <- tmp$customize.ring1
customize.ring2 <- tmp$customize.ring2
customize.ring3 <- tmp$customize.ring3
}
}
# LSC covert
if (lsc.converse){
if (nrow(ir_table) == 1){
warning("Didn't get IR region from thid species. ",
"It's impossible to convert LSC region.")
gene_table <- gene_table
} else{
tmp <- convert_region(ir_table = ir_table, l = l, region = "LSC",
genome = genome, gene_table = gene_table,
customize.ring1 = customize.ring1,
customize.ring2 = customize.ring2,
customize.ring3 = customize.ring3)
gene_table <- tmp$gene_table
gc_count <- tmp$gc_count
gc_total <- tmp$gc_total
customize.ring1 <- tmp$customize.ring1
customize.ring2 <- tmp$customize.ring2
customize.ring3 <- tmp$customize.ring3
}
}
# IRA covert
if (ira.converse){
if (nrow(ir_table) == 1){
warning("Didn't get IR region from thid species. ",
"It's impossible to convert IRA region.")
} else{
tmp <- convert_region(ir_table = ir_table, l = l, region = "IRA",
genome = genome, gene_table = gene_table,
customize.ring1 = customize.ring1,
customize.ring2 = customize.ring2,
customize.ring3 = customize.ring3,
indel_table = indel_table)
gene_table <- tmp$gene_table
gc_count <- tmp$gc_count
gc_total <- tmp$gc_total
customize.ring1 <- tmp$customize.ring1
customize.ring2 <- tmp$customize.ring2
customize.ring3 <- tmp$customize.ring3
indel_table <- tmp$indel_table
}
}
# IRB covert
if (irb.converse){
if (nrow(ir_table) == 1){
warning("Didn't get IR region from thid species. ",
"It's impossible to convert IRB region.")
gene_table <- gene_table
} else{
tmp <- convert_region(ir_table = ir_table, l = l, region = "IRB",
genome = genome, gene_table = gene_table,
customize.ring1 = customize.ring1,
customize.ring2 = customize.ring2,
customize.ring3 = customize.ring3,
indel_table = indel_table)
gene_table <- tmp$gene_table
gc_count <- tmp$gc_count
gc_total <- tmp$gc_total
customize.ring1 <- tmp$customize.ring1
customize.ring2 <- tmp$customize.ring2
customize.ring3 <- tmp$customize.ring3
indel_table <- tmp$indel_table
}
}
# colouring
color_table <- plastidGeneColor(psa.color = psa.color, psb.color = psb.color,
pet.color = pet.color, atp.color = atp.color,
ndh.color = ndh.color, rbc.color = rbc.color,
rpo.color = rpo.color, rsp.color = rsp.color,
rpl.color = rpl.color, clp_mat_inf.color = clp_mat_inf.color,
ycf.color = ycf.color, trn.color = trn.color,
rrn.color = rrn.color, other_gene.color = other_gene.color)
gene_table <- dplyr::left_join(gene_table, color_table, by = c("class" = "acronym"))
gene_table <- gene_table[which(gene_table$class %in% show.genes), ]
gene_table$gene <- gsub('(.{1,20})(\\s)', '\\1\n', gene_table$gene)
color_table <- color_table[which(color_table$acronym %in% unique(gene_table$class)), ]
# gene labels
gene_table$gene_label <- gene_table$gene
if (cu.bias){
gene_table$gene_label <- apply(gene_table[, c("gene_label", "cu_bias")], 1,
function(x){
if (is.na(x[2])){
return(x[1])
} else {
return(paste0(x[1], " (",
signif(as.numeric(x[2]), 2), ")"))
}
})
}
if (pseudo){
gene_table$gene_label <- apply(gene_table[, c("gene_label", "pseudo")], 1,
function(x){
if (x[2]){
return(paste0(x[1], " *"))
} else {
return(x[1])
}
})
}
genome <- genome
gene_table_r <- filter(gene_table, strand == "-") %>%
select(chr, start, end, gene, col, gc_count, gene_label) %>%
arrange(start)
gene_table_f <- filter(gene_table, strand == "+") %>%
select(chr, start, end, gene, col, gc_count, gene_label) %>%
arrange(start)
# 2. Set colors
ir_table$bg_col <- rep(NA, nrow(ir_table))
ir_table$bg_col[ir_table$name == "LSC"] <- lsc.color
ir_table$bg_col[ir_table$name == "SSC"] <- ssc.color
ir_table$bg_col[grepl("IR", ir_table$name)] <- ir.color
ir_table$bg_col[ir_table$name == "Genome"] <- lsc.color
ir_table$axis_col <- rep(NA, nrow(ir_table))
ir_table$axis_col[ir_table$name == "LSC"] <- gene_axis_lsc.color
ir_table$axis_col[ir_table$name == "SSC"] <- gene_axis_ssc.color
ir_table$axis_col[grepl("IR", ir_table$name)] <- gene_axis_ir.color
ir_table$axis_col[ir_table$name == "Genome"] <- gene_axis_lsc.color
# Automatically adjust colors
info.color <- CompColor(info.background)
ir_table$inf_col <- rep(NA, nrow(ir_table))
ir_table$inf_col[ir_table$name=="LSC"] <- CompColor(lsc.color)
ir_table$inf_col[ir_table$name=="SSC"] <- CompColor(ssc.color)
ir_table$inf_col[grepl("IR", ir_table$name)] <- CompColor(ir.color)
ir_table$inf_col[ir_table$name == "Genome"] <- CompColor(lsc.color)
if (!is.null(indel_table)){
indel_table$bg_color <- rep("White", #Transparent("#FFFFFF", 0.9),
nrow(indel_table))
}
ir_table$y <- rep(0.25, nrow(ir_table))
ir_table$y[grepl("IR", ir_table$name)] <- 0.75
# 3. Genome information
n_trn <- sum(grepl("trn", gene_table$gene, perl = TRUE))
n_rrn <- sum(grepl("rrn", gene_table$gene, perl = TRUE))
n_gene <- nrow(gene_table) - n_trn - n_rrn
# 4. Generate plot
# Initialize the plot device
if (save){
if (is.null(file.name)){
file <- sp_name
} else {
file <- file.name
}
if (!file.type %in% c("jpeg", "bmp", "png", "tiff", "pdf")){
warning("Can not save plot in ", file.type, " format. Avaliable formats
are 'jpeg', 'bmp', 'png', 'tiff',and 'pdf'.")
} else if (file.type == "pdf") {
grDevices::pdf(paste(file, file.type, sep="."), width=10.3, height=8.9)
} else {
do.call(file.type, args=list(filename=paste(file, file.type, sep="."),
width = 10.3, height = 8.9, units = "in",
res = 1000))
}
}
# Initialize the layout
ssc_deg<-ir_table$center[which(ir_table$name=="SSC")][1]
if ((length(ssc_deg) == 1) & !is.na(ssc_deg)){
ssc_deg <- 360 - round(((ssc_deg)/l) * 360)
rotate <- 270 - ssc_deg
} else {
rotate <- -90
}
# Track indexes
track_index <- 1:8
names(track_index) <- c("gene out 1", "gene out 2", "gene", "gene inner 1",
"gene inner 2", "arrow", "gc count", "ir region")
if (!is.null(customize.ring1)){
track_index[c("arrow", "gc count", "ir region")] <- track_index[c("arrow", "gc count", "ir region")] + 1
track_index["customize 1"] <- track_index["gene inner 2"] + 1
}
if (!is.null(customize.ring2)){
track_index[c("arrow", "gc count", "ir region")] <- track_index[c("arrow", "gc count", "ir region")] + 1
track_index["customize 2"] <- track_index["arrow"] - 1
}
if (!is.null(customize.ring3)){
track_index[c("arrow", "gc count", "ir region")] <- track_index[c("arrow", "gc count", "ir region")] + 1
track_index["customize 3"] <- track_index["arrow"] - 1
}
# Track heights
n_customize_ring <- sum(c(!is.null(customize.ring1), !is.null(customize.ring2), !is.null(customize.ring3)))
if (n_customize_ring == 1){
height[2] <- height[2]/1.5
} else if (n_customize_ring == 2) {
height[2] <- height[2]/2
height[3] <- height[3] * 0.7
} else if (n_customize_ring == 3) {
height[2] <- height[2]/2
height[3] <- height[3] * 0.7
height[1] <- height[1] *0.7
}
circlize::circos.clear()
circlize::circos.par(start.degree = rotate,
gap.after = 0,
track.margin = c(0, 0), cell.padding = c(0, 0, 0, 0))
circlize::circos.genomicInitialize(data=data.frame(chr="chr1", start=0, end=l,
stringsAsFactors = FALSE),
plotType = NULL)
# 1.2. gene label outside
if (nrow(gene_table_r) == 0 ) {
# Background
circlize::draw.sector(0, 360,
rou1 = 1,
rou2 = 0,
col = background, border = NA)
} else {
circlize::circos.genomicLabels(gene_table_r, labels = gene_table_r$gene_label,
side = "outside",cex = 0.5 * text.size,
connection_height = circlize::convert_height(3, "mm"),
labels_height = max(graphics::strwidth(gene_table_r$gene_label,
cex = 0.5 * text.size,
font = graphics::par("font")))
)
# Background
circlize::draw.sector(0, 360,
rou1 = circlize::get.cell.meta.data("cell.bottom.radius",
track.index = 2),
rou2 = 0,
col = background, border = NA)
}
# 3. Gene rectangles
if (nrow(gene_table_r) != 0 & nrow(gene_table_f) != 0) {
circlize::circos.genomicTrack(gene_table_f, factors = as.factor("chr1"),
ylim = c(-1, 1), bg.border = NA,
track.height = height[1], #circlize::convert_height(10, "mm"),
panel.fun = function(region, value, ...) {
if(gc.per.gene){
# genes lay on forward chain
circlize::circos.genomicRect(
gene_table_f[, c("start", "end")],
value = gene_table_f$gene,
ybottom = - gene_table_f$gc_count,
ytop = 0,
col = gene_table_f$col,
border = NA)
circlize::circos.genomicRect(
gene_table_f[, c("start", "end")],
value = gene_table_f$gene,
ybottom = -1,
ytop = 0,
lwd = 0.5,
col = Transparent(gene_table_f$col, 0.7),
border = NA
)
# genes lay on reverse chain
circlize::circos.genomicRect(
gene_table_r[, c("start", "end")],
value = gene_table_r$gene,
ybottom = 0,
ytop = gene_table_r$gc_count,
col = gene_table_r$col,
border = NA)
circlize::circos.genomicRect(
gene_table_r[, c("start", "end")],
value = gene_table_r$gene,
ybottom = 0,
ytop = 1,
lwd = 0.5,
col = Transparent(gene_table_r$col, 0.7),
border = NA
)
} else {
# genes lay on forward chain
circlize::circos.genomicRect(
gene_table_f[, c("start", "end")],
value = gene_table_f$gene,
ybottom = -1,
ytop = 0,
lwd = 0.5,
col = gene_table_f$col
)
circlize::circos.genomicRect(
gene_table_r[, c("start", "end")],
value = gene_table_r$gene,
ybottom = 0,
ytop = 1,
lwd = 0.5,
col = gene_table_r$col
)
}
# the bars in the middle indicate the IR, SSR, LSR
circlize::circos.rect(xleft = ir_table$start,
ybottom = - 0.05,
xright = ir_table$end,
ytop = 0.05,
col = ir_table$axis_col,
border = NA)
})
} else if (nrow(gene_table_f) != 0){
circlize::circos.genomicTrack(gene_table_f, factors = as.factor("chr1"),
ylim = c(-1, 1), bg.border = NA,
track.height = height[1], #circlize::convert_height(10, "mm"),
panel.fun = function(region, value, ...) {
if (gc.per.gene){
# genes lay on forward chain
circlize::circos.genomicRect(
gene_table_f[, c("start", "end")],
value = gene_table_f$gene,
ybottom = - gene_table_f$gc_count,
ytop = 0,
col = gene_table_f$col,
border = NA)
circlize::circos.genomicRect(
gene_table_f[, c("start", "end")],
value = gene_table_f$gene,
ybottom = -1,
ytop = 0,
lwd = 0.5,
col = Transparent(gene_table_f$col, 0.7),
border = NA
)
} else {
# genes lay on forward chain
circlize::circos.genomicRect(
gene_table_f[, c("start", "end")],
value = gene_table_f$gene,
ybottom = -1,
ytop = 0,
lwd = 0.5,
col = gene_table_f$col
)
}
# the bars in the middle indicate the IR, SSR, LSR
circlize::circos.rect(xleft = ir_table$start,
ybottom = - 0.05,
xright = ir_table$end,
ytop = 0.05,
col = ir_table$axis_col,
border = NA)
})
track_index[3:length(track_index)] <- track_index[3:length(track_index)] - 2
track_index <- track_index[-which(names(track_index) %in% c("gene out 1", "gene out 2"))]
} else if (nrow(gene_table_r) != 0) {
circlize::circos.genomicTrack(gene_table_r, factors = as.factor("chr1"),
ylim = c(-1, 1), bg.border = NA,
track.height = height[1], #circlize::convert_height(10, "mm"),
panel.fun = function(region, value, ...) {
if(gc.per.gene){
# genes lay on reverse chain
circlize::circos.genomicRect(
gene_table_r[, c("start", "end")],
value = gene_table_r$gene,
ybottom = 0,
ytop = gene_table_r$gc_count,
col = gene_table_r$col,
border = NA)
circlize::circos.genomicRect(
gene_table_r[, c("start", "end")],
value = gene_table_r$gene,
ybottom = 0,
ytop = 1,
lwd = 0.5,
col = Transparent(gene_table_r$col, 0.7),
border = NA
)
} else {
circlize::circos.genomicRect(
gene_table_r[, c("start", "end")],
value = gene_table_r$gene,
ybottom = 0,
ytop = 1,
lwd = 0.5,
col = gene_table_r$col
)
}
# the bars in the middle indicate the IR, SSR, LSR
circlize::circos.rect(xleft = ir_table$start,
ybottom = - 0.05,
xright = ir_table$end,
ytop = 0.05,
col = ir_table$axis_col,
border = NA)
})
track_index[(which(names(track_index) == "gene inner 2") + 1):length(track_index)] <-
track_index[(which(names(track_index) == "gene inner 2") + 1):length(track_index)] - 2
track_index <- track_index[-which(names(track_index) %in% c("gene inner 1", "gene inner 2"))]
} else {
track_index <- track_index[6:length(track_index)] -5
}
# 4.5. gene labels inside
if (nrow(gene_table_f) != 0){
circlize::circos.genomicLabels(gene_table_f, labels = gene_table_f$gene_label,
side = "inside", cex = 0.5 * text.size,
connection_height = circlize::convert_height(3, "mm"),
labels_height = max(graphics::strwidth(gene_table_f$gene_label,
cex = 0.5 * text.size,
font = graphics::par("font"))))
}
# customized ring 1
if (!is.null(customize.ring1)) {
customize.ring1 <- customize.ring1[order(customize.ring1$position), ]
customize.ring1$chr <- rep("chr1", nrow(customize.ring1))
ylim <- c(min(customize.ring1$value), max(customize.ring1$value))
bar_color <- ifelse(customize.ring1$value > 0, "orangered3", "steelblue")
if (min(customize.ring1$value) > 0) {
ylim[1] <- 0
bar_color <- rep(customize.ring1.color, nrow(customize.ring1))
}
if (max(customize.ring1$value) < 0){
ylim[2] <- 0
bar_color <- rep(customize.ring1.color, nrow(customize.ring1))
}
style <- switch(customize.ring1.type,
"line" = list(type = "l", arear = FALSE, col = customize.ring1.color, lwd = 0.5),
"line + filling" = list(type = "l", arear = TRUE, col = customize.ring1.color, lwd = 0.01),
"line + dot" = list(type = "o", arear = FALSE, col = bar_color, lwd = 0.5),
"line + dot + filling" = list(type = "o", arear = TRUE, col = bar_color, lwd = 0.01),
"step line" = list(type = "s", arear = FALSE, col = customize.ring1.color, lwd = 0.5),
"step line + filling" = list(type = "s", arear = TRUE, col = customize.ring1.color, lwd = 0.01),
"vertical line" = list(type = "h", arear = FALSE, col = bar_color, lwd = 0.5))
circlize::circos.track(factors =as.factor(customize.ring1$chr),
x = customize.ring1$position,
y = customize.ring1$value, ylim = ylim,
track.height = height[2]/2,
track.margin = c(0, circlize::convert_height(1, "mm")),
bg.border = NA, bg.col = gc.background,
panel.fun = function(x, y) {
circlize::circos.lines(x, y, type = style$type,
area = style$area,
baseline = 0,
pch = 1,
cex = 0.2 * text.size,
col = style$col,
lwd = style$lwd)
circlize::circos.segments(x0=0, x1=l, y0=0,
y1=0, lwd=0.5,
lty="16",
col=darken(customize.ring1.color, 0.7))
})
}
# customized ring 2
if (!is.null(customize.ring2)) {
customize.ring2 <- customize.ring2[order(customize.ring2$position), ]
customize.ring2$chr <- rep("chr1", nrow(customize.ring2))
ylim <- c(min(customize.ring2$value), max(customize.ring2$value))
bar_color <- ifelse(customize.ring2$value > 0, "orangered3", "steelblue")
if (min(customize.ring2$value) > 0) {
ylim[1] <- 0
bar_color <- rep(customize.ring2.color, nrow(customize.ring2))
}
if (max(customize.ring2$value) < 0){
ylim[2] <- 0
bar_color <- rep(customize.ring2.color, nrow(customize.ring2))
}
style <- switch(customize.ring2.type,
"line" = list(type = "l", arear = FALSE, col = customize.ring2.color, lwd = 0.5),
"line + filling" = list(type = "l", arear = TRUE, col = customize.ring2.color, lwd = 0.01),
"line + dot" = list(type = "o", arear = FALSE, col = bar_color, lwd = 0.5),
"line + dot + filling" = list(type = "o", arear = TRUE, col = bar_color, lwd = 0.01),
"step line" = list(type = "s", arear = FALSE, col = customize.ring2.color, lwd = 0.5),
"step line + filling" = list(type = "s", arear = TRUE, col = customize.ring2.color, lwd = 0.01),
"vertical line" = list(type = "h", arear = FALSE, col = bar_color, lwd = 0.5))
circlize::circos.track(factors =as.factor(customize.ring2$chr),
x = customize.ring2$position,
y = customize.ring2$value, ylim = ylim,
track.height = height[2]/2,
track.margin = c(0, circlize::convert_height(1, "mm")),
bg.border = NA, bg.col = gc.background,
panel.fun = function(x, y) {
circlize::circos.lines(x, y, type = style$type,
area = style$area,
baseline = 0,
pch = 1,
cex = 0.2 * text.size,
col = style$col,
lwd = style$lwd)
circlize::circos.segments(x0=0, x1=l, y0=0,
y1=0, lwd=0.5,
lty="16",
col=darken(customize.ring2.color, 0.7))
})
}
# customize.ring 3
if (!is.null(customize.ring3)) {
customize.ring3$chr <- rep("chr1", nrow(customize.ring3))
ylim <- c(min(customize.ring3$value), max(customize.ring3$value))
if (min(customize.ring3$value) > 0) {
ylim[1] <- 0
}
if (max(customize.ring3$value) < 0){
ylim[2] <- 0
}
circlize::circos.track(factors =as.factor(customize.ring3$chr),
track.margin = c(0, circlize::convert_height(1, "mm")),
track.height = height[2]/2,
ylim = ylim,
bg.border = NA, bg.col = gc.background,
panel.fun = function(x, y) {
if (min(customize.ring3$value) < 0 &
max(customize.ring3$value) > 0) {
df_large <- customize.ring3[which(customize.ring3$value > 0), ]
df_less <- customize.ring3[which(customize.ring3$value < 0), ]
circlize::circos.rect(xleft = df_large$start,
xright = df_large$end,
ybottom = 0,
ytop = df_large$value,
col = "orangered3",
border = NA)
circlize::circos.rect(xleft = df_less$start,
xright = df_less$end,
ybottom = df_less$value,
ytop = 0,
col = "steelblue",
border = NA)
} else if (min(customize.ring3$value) < 0) {
circlize::circos.rect(xleft = customize.ring3$start,
xright = customize.ring3$end,
ybottom = customize.ring3$value,
ytop = 0,
col = customize.ring3.color,
border = NA)
} else {
circlize::circos.rect(xleft = customize.ring3$start,
xright = customize.ring3$end,
ybottom = 0,
ytop = customize.ring3$value,
col = customize.ring3.color,
border = NA)
}
})
}
# 6. Arrow outside gnome axis
circlize::circos.track(ylim = c(0, 1), track.height = 0.06, bg.border = NA,
track.margin = c(0, circlize::convert_height(1, "mm")),
panel.fun = function(x, y) {
circlize::circos.arrow(x1 = l - l %/% 40,
x2 = l,
y = 0.85,
arrow.position="start",
col = MildColor(background),
border = NA,
arrow.head.length = circlize::convert_x(3, "mm"),
width = 0.15)
# circlize::circos.genomicAxis(h = "bottom")
})
# 7. GC content
circlize::circos.track(factors =as.factor(gc_count$chr),
x=gc_count$position,
y = 1- gc_count$gc_count,
ylim = c(0, 1),
track.height = height[2],
bg.border = NA, bg.col = gc.background,
panel.fun = function(x, y) {
circlize::circos.lines(x, y, type = "s", area = TRUE,
baseline = "top",
col = gc.color, lwd = 0.00001)
circlize::circos.arrow(x1= 0,
x2= l %/% 35,
y= 0.9,
col=MildColor(gc.color),
border=NA,
arrow.head.length=circlize::convert_x(3, "mm"),
width=0.15 * 0.05 / height[2])
circlize::circos.segments(x0=0, x1=l, y0=0.25,
y1=0.25, lwd=0.5,
lty="16", col=MildColor(gc.background))
circlize::circos.segments(x0=0, x1=l, y0=0.5,
y1=0.5, lwd=0.6,
lty="11", col=MildColor(gc.background))
circlize::circos.segments(x0=0, x1=l, y0=0.75,
y1=0.75, lwd=0.5,
lty="16", col=MildColor(gc.color))
circlize::circos.genomicAxis(h = "top", labels.cex = 0.4 * text.size)
})
#circos.yaxis(at=c(0.25,0.5,0.75),labels.cex=0.25 * text.size,lwd=0,tick.length=0,
# labels.col="grey40",col="#FFFFFF")
# 8. inner ring for IR, SSR, LSR
circlize::circos.genomicTrack(ir_table, bg.border = NA,
ylim = c(0, 1),track.height = height[3],
panel.fun = function(region, value, ...) {
if (ir.gc){
circlize::circos.rect(xleft = ir_table$start,
ybottom = 0,#circlize::convert_y(0, "mm"),
xright = ir_table$end,
ytop = 1,#circlize::convert_y(3, "mm"),
col = Transparent(ir_table$bg_col, 0.7),
border = NA)
circlize::circos.rect(xleft = ir_table$start,
ybottom =0,#circlize::convert_y(0, "mm"),
xright = ir_table$end,
ytop = ir_table$gc_count,#circlize::convert_y(3, "mm"),
col = ir_table$bg_col,
border = NA)
} else {
circlize::circos.rect(xleft = ir_table$start,
ybottom = 0,#circlize::convert_y(0, "mm"),
xright = ir_table$end,
ytop = 1,#circlize::convert_y(3, "mm"),
col = ir_table$bg_col,
border = NA)
}
if (!is.null(indel_table) & show.indel){
circlize::circos.rect(xleft = indel_table$position - 10,
ybottom = 0,#circlize::convert_y(0, "mm"),
xright = indel_table$position + 10,
ytop = 1,#circlize::convert_y(3, "mm"),
col = indel_table$bg_col,
border = NA)
circlize::circos.points(x = indel_table$position,
y = 0.05,
col = indel_table$col,
pch = 19,
cex = 0.5 * text.size)
# circlize::circos.text(x = indel_table$position,
# y = 0.5,
# col = indel_table$col,
# labels = indel_table$string,
# cex = 0.2 * text.size)
}
circlize::circos.text(x = ir_table$center,
y = ir_table$y, cex = 0.7 * text.size,
labels = ir_table$text,
facing = "bending.inside",
col = ir_table$inf_col,
niceFacing = TRUE)
})
# Specie's information inthe central
circlize::draw.sector(0, 360,
rou1 = circlize::get.cell.meta.data("cell.bottom.radius",
track.index = track_index['ir region']),
rou2 = 0,
col = info.background, border = NA)
if (organelle_type){
graphics::text(0,0.10, sp_name, font=4, cex = 0.8 * text.size, col = info.color)
graphics::text(0,0.05, "Plastid Genome", font=4, cex = 0.8 * text.size, col = info.color)
} else{
graphics::text(0,0.05, sp_name, font=4, cex = 0.8 * text.size, col = info.color)
}
ns <- c(paste(n_gene, "genes"), paste(n_rrn, "rRNAs"), paste(n_trn, "tRNAs"))[which(c(gene.no, rrn.no, trn.no))]
text <- c(paste(ns, collapse = "; "), paste(prettyNum(l, big.mark = ","), "bp", " "),
paste("GC:",round(gc_total, 2)*100,"%"))
j <- 1
if (cu.bias){
cex = 0.7 * text.size
} else {
cex = 0.8 * text.size
}
for (i in which(c((length(ns) != 0), genome.length, total.gc))) {
graphics::text(0,0.05 - 0.05 * j, text[i],
font=4, cex = cex, col = info.color)
j <- j + 1
}
# Highlight IR
if (shadow & (nrow(ir_table) >= 4)){
pos_s=circlize::circlize(ir_table$start[which(grepl("IR",
ir_table$name))],
1, sector.index = "chr1", track.index = 1)
pos_e=circlize::circlize(ir_table$end[which(grepl("IR",
ir_table$name))],
1, sector.index = "chr1", track.index = 1)
for (i in 1:nrow(pos_s)) {
if ('gene' %in% names(track_index)){
circlize::draw.sector(pos_s[i, "theta"], pos_e[i, "theta"],
#start.degree = 0, end.degree = 360 - 10,
rou1 = circlize::get.cell.meta.data("cell.top.radius",
track.index = track_index['gene']),
rou2 = circlize::get.cell.meta.data("cell.bottom.radius",
track.index = track_index['gc count']),
col = shadow.color, border = NA)
} else {
circlize::draw.sector(pos_s[i, "theta"], pos_e[i, "theta"],
#start.degree = 0, end.degree = 360 - 10,
rou1 = circlize::get.cell.meta.data("cell.top.radius",
track.index = track_index['arrow']),
rou2 = circlize::get.cell.meta.data("cell.bottom.radius",
track.index = track_index['gc count']),
col = shadow.color, border = NA)
}
}
if (!is.null(indel_table) & show.indel) {
pos_s=circlize::circlize(indel_table$position - 10,
1, sector.index = "chr1", track.index = 1)
pos_e=circlize::circlize(indel_table$position + 10,
1, sector.index = "chr1", track.index = 1)
for (i in 1:nrow(pos_s)) {
if ('gene' %in% names(track_index)){
circlize::draw.sector(pos_s[i, "theta"], pos_e[i, "theta"],
#start.degree = 0, end.degree = 360 - 10,
rou1 = circlize::get.cell.meta.data("cell.top.radius",
track.index = track_index['gene']),
rou2 = circlize::get.cell.meta.data("cell.bottom.radius",
track.index = track_index['gc count']),
col = Transparent("#FFFFFF", 0.5),
border = NA)
} else {
circlize::draw.sector(pos_s[i, "theta"], pos_e[i, "theta"],
#start.degree = 0, end.degree = 360 - 10,
rou1 = circlize::get.cell.meta.data("cell.top.radius",
track.index = track_index['arrow']),
rou2 = circlize::get.cell.meta.data("cell.bottom.radius",
track.index = track_index['gc count']),
col = Transparent("#FFFFFF", 0.9),
border = NA)
}
}
}
}
# Add legend
if (legend){
legend(x = -1.27, y = -0.6, legend = color_table$label, cex = 0.6 * text.size,
fill = color_table$col, box.col = "white")
}
if (save){
grDevices::dev.off()
}
}
#' Generate Mitochondrion Genome Plot
#'
#' @param plot.tables a list contains information of IR region, gene, and gc
#' count of the genome. It can be generated by function \code{\link{PlotTab}}.
#' @param save A logical value. If it is \code{TRUE}, the plot will be saved in
#' work directory.
#' @param file.type A charactor. It indicates the format of the file in which
#' the plot will be saved. Options are: "pdf", "png", "jpeg","bmp", "tiff".
#' By default, it is set as "pdf".
#' @param file.name A charactor. It indicates the name of the file in which
#' the plot will be saved. By default, it is set as the specie's name.
#' @param legend A logical value. If it is \code{TRUE}, the legend for gene
#' colors will be shown.
#' @param background An R color object. It indicates the color for the background of
#' entire plot area.
#' @param height A vector of numeric value. The elements of it indicat the
#' height of gene layer, GC count layer and IR region layer, respectively.
#' Default setting is "0.1, 0.2, 0.07". The taltal circle plot region always
#' has a radius of 1, so a height of 0.1 means 10\% of the circle radius.
#' @param gc.color An R color object. It indicates the color for the lines in gc count
#' plot.
#' @param gc.background An R color object. It indicates the color for the background
#' of gc count plot area.
#' @param info.background An R color object. It indicates the color for the background
#' of central area where species' information was shown.
#' @param gc.per.gene A logical value. If it is \code{TRUE}, the gc content for
#' each gene will be plot on the gene layer with deep colors.
#' @param pseudo A logical value. If it is \code{TRUE}, the pseudo genes (if
#' there is some in the species) will be marked with a "*" at the end of the
#' labels.
#' @param cu.bias A logical value. If it is \code{TRUE}, the condon usage bias
#' metric for each gene will be shown in the labels. The metric “Measure
#' Independent of Length and Composition (MILC)” was used for evaluate the bias.
#' You can find more details in the reference paper.
#' @param text.size A numeric value. It indicates the size of all texts in the
#' plot. For exmple, \code{text.size = 1.5} means all the texts in the plot will
#' be enlarged as 1.5 times of their original size.
#' @param gc.per.gene A logical value. If it is \code{TRUE}, the GC content of
#' each gene will be show by a darker part in gene rectangles.
#' @param genome.length A logical value. If it is \code{TRUE}, the length of
#' genome will be shown in the center of the plot.
#' @param total.gc A logical value. If it is \code{TRUE}, the GC content of
#' whole genome will be shown in the center of the plot.
#' @param gene.no A logical value. If it is \code{TRUE}, the number of genes in
#' the genome will be shown in the center of the plot.
#' @param rrn.no A logical value. If it is \code{TRUE}, the number of rRNAs in
#' the genome will be shown in the center of the plot.
#' @param trn.no A logical value. If it is \code{TRUE}, the number of tRNAs in
#' the genome will be shown in the center of the plot.
#' @param organelle A logical value. If it is \code{TRUE}, the organelle type of
#' the genome will be shown in the center of the plot.
#' @param nad.color An R color string. It indicates the color for genes of
#' complex I (NADH dehydrogenase)
#' @param sdh.color An R color string. It indicates the color for genes of
#' complex II (succinate dehydrogenase)
#' @param cob.color An R color string. It indicates the color for genes of
#' complex III (ubichinol cytochrome reductase)
#' @param cox.color An R color string. It indicates the color for genes of
#' complex IV (cytochrome c oxidase)
#' @param atp.color An R color string. It indicates the color for genes of
#' ATP synthase
#' @param ccmF.color An R color string. It indicates the color for genes of
#' cytochrome c biogenesi
#' @param mmt.color An R color string. It indicates the color for genes of
#' transport membrane protein
#' @param rps.color An R color string. It indicates the color for genes of
#' ribosomal proteins (SSU)
#' @param rpl.color An R color string. It indicates the color for genes of
#' ribosomal proteins (LSU)
#' @param mat.color An R color string. It indicates the color for genes of
#' maturases
#' @param orf.color An R color string. It indicates the color for ORFs
#' @param trn.color An R color object. It indicates the color for genes of
#' transfer RNA
#' @param rrn.color An R color object. It indicates the color for genes of
#' ribosomal RNA
#' @param other_gene.color An R color object. It indicates the color for other
#' genes
#' @param show.gene A vector of characters. It indicates which classes of genes
#' will be shown on the plot. Avaliable values are "nad|nd","sdh","cob",
#' "cox|cytb","atp", "ccmF","mtt","rps","rpl", "mat","orf","trn","rrn", "OTHER"
#' @param customize.ring1 A data frame. It must contain 2 columns:
#' \itemize{
#' \item \strong{position}: 1-base genomic coordinate for the features.
#' \item \strong{value}: the values for the features.
#' }
#'
#' @param customize.ring1.type A character. It indicate the plot type in
#' customize.ring1. Available values are "line", "line + filling", "line + dot",
#' "line + dot + filling", "step line", "step line + filling", "vertical line"
#' @param customize.ring1.color An R color object. It indicates the color for
#' the plots in customized ring 1.
#' @param customize.ring2 A data frame. It must contain 2 columns:
#' \itemize{
#' \item \strong{position}: 1-base genome coordinate for the features.
#' \item \strong{value}: the values for the features.
#' }
#' @param customize.ring2.type A character. It indicate the plot type in
#' customize.ring2. Available values are "line", "line + filling", "line + dot",
#' "line + dot + filling", "step line", "step line + filling", "vertical line"
#' @param customize.ring2.color An R color object. It indicates the color for
#' the plots in customized ring 2.
#' @param customize.ring3 A data frame. It must contain 2 columns:
#' \itemize{
#' \item \strong{start}: 1-base genome coordinate for the start point of the
#' features.
#' \item \strong{end}: 1-base genome coordinate for the end point of the
#' features.
#' \item \strong{value}: the values for the features.
#' }
#' @param customize.ring3.color An R color object. It indicates the color for
#' the plots in customized ring 3.
#' @return A plot for chloroplast genome.
#'
#' @references Supek, Fran, and Kristian Vlahovicek. “Comparison of codon usage
#' measures and their applicability in prediction of microbial gene
#' expressivity.” BMC bioinformatics vol. 6 182. 19 Jul.
#' 2005, doi:10.1186/1471-2105-6-182
#'
#' @importFrom magrittr %>%
#' @import dplyr
#' @export
PlotMitGenome <- function(plot.tables, save = TRUE, file.type = "pdf",
text.size = 1, height = c(0.1, 0.2, 0.07),
file.name = NULL,
gc.per.gene = TRUE, pseudo = TRUE, legend = TRUE,
genome.length = TRUE, total.gc = TRUE,
gene.no = TRUE, rrn.no = TRUE, trn.no = TRUE,
organelle_type = TRUE,
background = "grey90", gc.color = "grey30",
gc.background = "grey70", info.background = "black",
nad.color = "#2A6332", sdh.color = "#4C8805",
cob.color = "#7F992C", cox.color = "#FEEE50",
atp.color = "#9FBB3D", ccmF.color = "#4D9E3F",
mmt.color = "#AE2D29", rps.color = "#D6AD7C",
rpl.color = "#9C7A4B", mat.color = "#D9662D",
orf.color = "#71B8A9", trn.color = "#172C7F",
rrn.color = "#D1382A", other_gene.color = "#7D7D7D",
show.genes = c("nad|nd","sdh","cob","cox|cytb","atp",
"ccmF","mtt","rps","rpl",
"mat","orf","trn","rrn", "OTHER"),
cu.bias = TRUE, customize.ring1 = NULL,
customize.ring1.color = "grey30",
customize.ring1.type = "line",
customize.ring2.color = "grey30",
customize.ring2 = NULL,
customize.ring2.type = "line",
customize.ring3 = NULL,
customize.ring3.color = "grey30"){
# Unpack plot.table
genome <- plot.tables$genome
gc_count <- plot.tables$gc_count
gc_total <- plot.tables$gc_total
sp_name <- plot.tables$sp_name
gene_table <- plot.tables$gene_table
genome <- plot.tables$genome
l <- Biostrings::nchar(genome)
# 1. Modify gene table
# colouring
color_table <- mitGeneColor(nad.color = nad.color, sdh.color = sdh.color,
cob.color = cob.color, cox.color = cox.color,
atp.color = atp.color, ccmF.color = ccmF.color,
mmt.color = mmt.color, rps.color = rps.color,
rpl.color = mmt.color, mat.color = mat.color,
orf.color = orf.color, trn.color = trn.color,
rrn.color = rrn.color,
other_gene.color = other_gene.color)
gene_table <- dplyr::left_join(gene_table, color_table, by = c("class" = "acronym"))
gene_table <- gene_table[which(gene_table$class %in% show.genes), ]
color_table <- color_table[which(color_table$acronym %in% unique(gene_table$class)), ]
# gene labels
gene_table$gene_label <- gene_table$gene
if (cu.bias){
gene_table$gene_label <- apply(gene_table[, c("gene_label", "cu_bias")], 1,
function(x){
if (is.na(x[2])){
return(x[1])
} else {
return(paste0(x[1], " (",
signif(as.numeric(x[2]), 2), ")"))
}
})
}
if (pseudo){
gene_table$gene_label <- apply(gene_table[, c("gene_label", "pseudo")], 1,
function(x){
if (x[2]){
return(paste0(x[1], " *"))
} else {
return(x[1])
}
})
}
genome <- genome
gene_table_r <- filter(gene_table, strand == "-") %>%
select(chr, start, end, gene, col, gc_count, gene_label) %>%
arrange(start)
gene_table_f <- filter(gene_table, strand == "+") %>%
select(chr, start, end, gene, col, gc_count, gene_label) %>%
arrange(start)
# 2. Set colors
# Automatically adjust colors
info.color <- CompColor(info.background)
# 3. Genome information
n_trn <- sum(grepl("trn", gene_table$gene, perl = TRUE))
n_rrn <- sum(grepl("rrn", gene_table$gene, perl = TRUE))
n_gene <- nrow(gene_table) - n_trn - n_rrn
# 4. Generate plot
# Initialize the plot device
if (save){
if (is.null(file.name)){
file <- sp_name
} else {
file <- file.name
}
if (!file.type %in% c("jpeg", "bmp", "png", "tiff", "pdf")){
warning("Can not save plot in ", file.type, " format. Avaliable formats
are 'jpeg', 'bmp', 'png', 'tiff',and 'pdf'.")
} else if (file.type == "pdf") {
grDevices::pdf(paste(file, file.type, sep="."), width=10.3, height=8.9)
} else {
do.call(file.type, args=list(filename=paste(file, file.type, sep="."),
width = 10.3, height = 8.9, units = "in",
res = 1000))
}
}
# Initialize the layout
# Track indexes
track_index <- 1:7
names(track_index) <- c("gene out 1", "gene out 2", "gene", "gene inner 1",
"gene inner 2", "arrow", "gc count")
if (!is.null(customize.ring1)){
track_index[c("arrow", "gc count")] <- track_index[c("arrow", "gc count")] + 1
track_index["customize 1"] <- track_index["gene inner 2"] + 1
}
if (!is.null(customize.ring2)){
track_index[c("arrow", "gc count")] <- track_index[c("arrow", "gc count")] + 1
track_index["customize 2"] <- track_index["arrow"] - 1
}
if (!is.null(customize.ring3)){
track_index[c("arrow", "gc count")] <- track_index[c("arrow", "gc count")] + 1
track_index["customize 3"] <- track_index["arrow"] - 1
}
# Track heights
n_customize_ring <- sum(c(!is.null(customize.ring1), !is.null(customize.ring2), !is.null(customize.ring3)))
if (n_customize_ring == 1){
height[2] <- height[2]/1.5
} else if (n_customize_ring == 2) {
height[2] <- height[2]/2
height[3] <- height[3] * 0.7
} else if (n_customize_ring == 3) {
height[2] <- height[2]/2
height[3] <- height[3] * 0.7
height[1] <- height[1] *0.7
}
circlize::circos.clear()
circlize::circos.par(start.degree = 0, gap.after = 0,
track.margin = c(0, 0), cell.padding = c(0, 0, 0, 0))
circlize::circos.genomicInitialize(data=data.frame(chr="chr1", start=0, end=l,
stringsAsFactors = FALSE),
plotType = NULL)
# 1.2. gene label outside
if (nrow(gene_table_r) == 0 ) {
# Background
circlize::draw.sector(0, 360,
rou1 = 1,
rou2 = 0,
col = background, border = NA)
} else {
circlize::circos.genomicLabels(gene_table_r, labels = gene_table_r$gene_label,
side = "outside",cex = 0.5 * text.size,
connection_height = circlize::convert_height(3, "mm"),
labels_height = max(graphics::strwidth(gene_table_r$gene_label,
cex = 0.5 * text.size,
font = graphics::par("font")))
)
# Background
circlize::draw.sector(0, 360,
rou1 = circlize::get.cell.meta.data("cell.bottom.radius",
track.index = 2),
rou2 = 0,
col = background, border = NA)
}
# 3. Gene rectangles
if (nrow(gene_table_r) != 0 & nrow(gene_table_f) != 0) {
circlize::circos.genomicTrack(gene_table_f, factors = as.factor("chr1"),
ylim = c(-1, 1), bg.border = NA,
track.height = height[1], #circlize::convert_height(10, "mm"),
panel.fun = function(region, value, ...) {
if(gc.per.gene){
# genes lay on forward chain
circlize::circos.genomicRect(
gene_table_f[, c("start", "end")],
value = gene_table_f$gene,
ybottom = - gene_table_f$gc_count,
ytop = 0,
col = gene_table_f$col,
border = NA)
circlize::circos.genomicRect(
gene_table_f[, c("start", "end")],
value = gene_table_f$gene,
ybottom = -1,
ytop = 0,
lwd = 0.5,
col = Transparent(gene_table_f$col, 0.7),
border = NA
)
# genes lay on reverse chain
circlize::circos.genomicRect(
gene_table_r[, c("start", "end")],
value = gene_table_r$gene,
ybottom = 0,
ytop = gene_table_r$gc_count,
col = gene_table_r$col,
border = NA)
circlize::circos.genomicRect(
gene_table_r[, c("start", "end")],
value = gene_table_r$gene,
ybottom = 0,
ytop = 1,
lwd = 0.5,
col = Transparent(gene_table_r$col, 0.7),
border = NA
)
} else {
# genes lay on forward chain
circlize::circos.genomicRect(
gene_table_f[, c("start", "end")],
value = gene_table_f$gene,
ybottom = -1,
ytop = 0,
lwd = 0.5,
col = gene_table_f$col
)
circlize::circos.genomicRect(
gene_table_r[, c("start", "end")],
value = gene_table_r$gene,
ybottom = 0,
ytop = 1,
lwd = 0.5,
col = gene_table_r$col
)
}
# the x axis
circlize::circos.rect(xleft = 0,
ybottom = - 0.05,
xright = l,
ytop = 0.05,
col = "black",
border = NA)
})
} else if (nrow(gene_table_f) != 0){
circlize::circos.genomicTrack(gene_table_f, factors = as.factor("chr1"),
ylim = c(-1, 1), bg.border = NA,
track.height = height[1], #circlize::convert_height(10, "mm"),
panel.fun = function(region, value, ...) {
if (gc.per.gene){
# genes lay on forward chain
circlize::circos.genomicRect(
gene_table_f[, c("start", "end")],
value = gene_table_f$gene,
ybottom = - gene_table_f$gc_count,
ytop = 0,
col = gene_table_f$col,
border = NA)
circlize::circos.genomicRect(
gene_table_f[, c("start", "end")],
value = gene_table_f$gene,
ybottom = -1,
ytop = 0,
lwd = 0.5,
col = Transparent(gene_table_f$col, 0.7),
border = NA
)
} else {
# genes lay on forward chain
circlize::circos.genomicRect(
gene_table_f[, c("start", "end")],
value = gene_table_f$gene,
ybottom = -1,
ytop = 0,
lwd = 0.5,
col = gene_table_f$col
)
}
# the x axis
circlize::circos.rect(xleft = 0,
ybottom = - 0.05,
xright = l,
ytop = 0.05,
col = "black",
border = NA)
})
track_index[3:length(track_index)] <- track_index[3:length(track_index)] - 2
track_index <- track_index[-which(names(track_index) %in% c("gene out 1", "gene out 2"))]
} else if (nrow(gene_table_r) != 0) {
circlize::circos.genomicTrack(gene_table_r, factors = as.factor("chr1"),
ylim = c(-1, 1), bg.border = NA,
track.height = height[1], #circlize::convert_height(10, "mm"),
panel.fun = function(region, value, ...) {
if(gc.per.gene){
# genes lay on reverse chain
circlize::circos.genomicRect(
gene_table_r[, c("start", "end")],
value = gene_table_r$gene,
ybottom = 0,
ytop = gene_table_r$gc_count,
col = gene_table_r$col,
border = NA)
circlize::circos.genomicRect(
gene_table_r[, c("start", "end")],
value = gene_table_r$gene,
ybottom = 0,
ytop = 1,
lwd = 0.5,
col = Transparent(gene_table_r$col, 0.7),
border = NA
)
} else {
circlize::circos.genomicRect(
gene_table_r[, c("start", "end")],
value = gene_table_r$gene,
ybottom = 0,
ytop = 1,
lwd = 0.5,
col = gene_table_r$col
)
}
# the x axis
circlize::circos.rect(xleft = 0,
ybottom = - 0.05,
xright = l,
ytop = 0.05,
col = "black",
border = NA)
})
track_index[(which(names(track_index) == "gene inner 2") + 1):length(track_index)] <-
track_index[(which(names(track_index) == "gene inner 2") + 1):length(track_index)] - 2
track_index <- track_index[-which(names(track_index) %in% c("gene inner 1", "gene inner 2"))]
} else {
track_index <- track_index[6:length(track_index)] -5
}
# 4.5. gene labels inside
if (nrow(gene_table_f) != 0){
circlize::circos.genomicLabels(gene_table_f, labels = gene_table_f$gene_label,
side = "inside", cex = 0.5 * text.size,
connection_height = circlize::convert_height(3, "mm"),
labels_height = max(graphics::strwidth(gene_table_f$gene_label,
cex = 0.5 * text.size,
font = graphics::par("font"))))
}
# customized ring 1
if (!is.null(customize.ring1)) {
customize.ring1 <- customize.ring1[order(customize.ring1$position), ]
customize.ring1$chr <- rep("chr1", nrow(customize.ring1))
ylim <- c(min(customize.ring1$value), max(customize.ring1$value))
bar_color <- ifelse(customize.ring1$value > 0, "orangered3", "steelblue")
if (min(customize.ring1$value) > 0) {
ylim[1] <- 0
bar_color <- rep(customize.ring1.color, nrow(customize.ring1))
}
if (max(customize.ring1$value) < 0){
ylim[2] <- 0
bar_color <- rep(customize.ring1.color, nrow(customize.ring1))
}
style <- switch(customize.ring1.type,
"line" = list(type = "l", arear = FALSE, col = customize.ring1.color, lwd = 0.5),
"line + filling" = list(type = "l", arear = TRUE, col = customize.ring1.color, lwd = 0.01),
"line + dot" = list(type = "o", arear = FALSE, col = bar_color, lwd = 0.5),
"line + dot + filling" = list(type = "o", arear = TRUE, col = bar_color, lwd = 0.01),
"step line" = list(type = "s", arear = FALSE, col = customize.ring1.color, lwd = 0.5),
"step line + filling" = list(type = "s", arear = TRUE, col = customize.ring1.color, lwd = 0.01),
"vertical line" = list(type = "h", arear = FALSE, col = bar_color, lwd = 0.5))
circlize::circos.track(factors =as.factor(customize.ring1$chr),
x = customize.ring1$position,
y = customize.ring1$value, ylim = ylim,
track.height = height[2]/2,
track.margin = c(0, circlize::convert_height(1, "mm")),
bg.border = NA, bg.col = gc.background,
panel.fun = function(x, y) {
circlize::circos.lines(x, y, type = style$type,
area = style$area,
baseline = 0,
pch = 1,
cex = 0.2 * text.size,
col = style$col,
lwd = style$lwd)
circlize::circos.segments(x0=0, x1=l, y0=0,
y1=0, lwd=0.5,
lty="16",
col=darken(customize.ring1.color, 0.7))
})
}
# customized ring 2
if (!is.null(customize.ring2)) {
customize.ring2 <- customize.ring2[order(customize.ring2$position), ]
customize.ring2$chr <- rep("chr1", nrow(customize.ring2))
ylim <- c(min(customize.ring2$value), max(customize.ring2$value))
bar_color <- ifelse(customize.ring2$value > 0, "orangered3", "steelblue")
if (min(customize.ring2$value) > 0) {
ylim[1] <- 0
bar_color <- rep(customize.ring2.color, nrow(customize.ring2))
}
if (max(customize.ring2$value) < 0){
ylim[2] <- 0
bar_color <- rep(customize.ring2.color, nrow(customize.ring2))
}
style <- switch(customize.ring2.type,
"line" = list(type = "l", arear = FALSE, col = customize.ring2.color, lwd = 0.5),
"line + filling" = list(type = "l", arear = TRUE, col = customize.ring2.color, lwd = 0.01),
"line + dot" = list(type = "o", arear = FALSE, col = bar_color, lwd = 0.5),
"line + dot + filling" = list(type = "o", arear = TRUE, col = bar_color, lwd = 0.01),
"step line" = list(type = "s", arear = FALSE, col = customize.ring2.color, lwd = 0.5),
"step line + filling" = list(type = "s", arear = TRUE, col = customize.ring2.color, lwd = 0.01),
"vertical line" = list(type = "h", arear = FALSE, col = bar_color, lwd = 0.5))
circlize::circos.track(factors =as.factor(customize.ring2$chr),
x = customize.ring2$position,
y = customize.ring2$value, ylim = ylim,
track.height = height[2]/2,
track.margin = c(0, circlize::convert_height(1, "mm")),
bg.border = NA, bg.col = gc.background,
panel.fun = function(x, y) {
circlize::circos.lines(x, y, type = style$type,
area = style$area,
baseline = 0,
pch = 1,
cex = 0.2 * text.size,
col = style$col,
lwd = style$lwd)
circlize::circos.segments(x0=0, x1=l, y0=0,
y1=0, lwd=0.5,
lty="16",
col=darken(customize.ring2.color, 0.7))
})
}
# customize.ring 3
if (!is.null(customize.ring3)) {
customize.ring3$chr <- rep("chr1", nrow(customize.ring3))
ylim <- c(min(customize.ring3$value), max(customize.ring3$value))
if (min(customize.ring3$value) > 0) {
ylim[1] <- 0
}
if (max(customize.ring3$value) < 0){
ylim[2] <- 0
}
circlize::circos.track(factors =as.factor(customize.ring3$chr),
track.margin = c(0, circlize::convert_height(1, "mm")),
track.height = height[2]/2,
ylim = ylim,
bg.border = NA, bg.col = gc.background,
panel.fun = function(x, y) {
if (min(customize.ring3$value) < 0 &
max(customize.ring3$value) > 0) {
df_large <- customize.ring3[which(customize.ring3$value > 0), ]
df_less <- customize.ring3[which(customize.ring3$value < 0), ]
circlize::circos.rect(xleft = df_large$start,
xright = df_large$end,
ybottom = 0,
ytop = df_large$value,
col = "orangered3",
border = NA)
circlize::circos.rect(xleft = df_less$start,
xright = df_less$end,
ybottom = df_less$value,
ytop = 0,
col = "steelblue",
border = NA)
} else if (min(customize.ring3$value) < 0) {
circlize::circos.rect(xleft = customize.ring3$start,
xright = customize.ring3$end,
ybottom = customize.ring3$value,
ytop = 0,
col = customize.ring3.color,
border = NA)
} else {
circlize::circos.rect(xleft = customize.ring3$start,
xright = customize.ring3$end,
ybottom = 0,
ytop = customize.ring3$value,
col = customize.ring3.color,
border = NA)
}
})
}
# 6. Arrow outside gnome axis
circlize::circos.track(ylim = c(0, 1), track.height = 0.06, bg.border = NA,
track.margin = c(0, circlize::convert_height(1, "mm")),
panel.fun = function(x, y) {
circlize::circos.arrow(x1 = l - l %/% 40,
x2 = l,
y = 0.85,
arrow.position="start",
col = MildColor(background),
border = NA,
arrow.head.length = circlize::convert_x(3, "mm"),
width = 0.15)
# circlize::circos.genomicAxis(h = "bottom")
})
# 7. GC content
circlize::circos.track(factors =as.factor(gc_count$chr),
x=gc_count$position,
y = 1- gc_count$gc_count,
ylim = c(0, 1),
track.height = height[2],
bg.border = NA, bg.col = gc.background,
panel.fun = function(x, y) {
circlize::circos.lines(x, y, type = "s", area = TRUE,
baseline = "top",
col = gc.color, lwd = 0.00001)
circlize::circos.arrow(x1= 0,
x2= l %/% 35,
y= 0.9,
col=MildColor(gc.color),
border=NA,
arrow.head.length=circlize::convert_x(3, "mm"),
width=0.15 * 0.05 / height[2])
circlize::circos.segments(x0=0, x1=l, y0=0.25,
y1=0.25, lwd=0.5,
lty="16", col=MildColor(gc.background))
circlize::circos.segments(x0=0, x1=l, y0=0.5,
y1=0.5, lwd=0.6,
lty="11", col=MildColor(gc.background))
circlize::circos.segments(x0=0, x1=l, y0=0.75,
y1=0.75, lwd=0.5,
lty="16", col=MildColor(gc.color))
circlize::circos.genomicAxis(h = "top", labels.cex = 0.4 * text.size)
})
#circos.yaxis(at=c(0.25,0.5,0.75),labels.cex=0.25 * text.size,lwd=0,tick.length=0,
# labels.col="grey40",col="#FFFFFF")
# Specie's information inthe central
circlize::draw.sector(0, 360,
rou1 = circlize::get.cell.meta.data("cell.bottom.radius",
track.index = track_index['gc count']),
rou2 = 0,
col = info.background, border = NA)
if (organelle_type){
graphics::text(0,0.10, sp_name, font=4, cex = 0.8 * text.size, col = info.color)
graphics::text(0,0.05, "Plastid Genome", font=4, cex = 0.8 * text.size, col = info.color)
} else{
graphics::text(0,0.05, sp_name, font=4, cex = 0.8 * text.size, col = info.color)
}
ns <- c(paste(n_gene, "genes"), paste(n_rrn, "rRNAs"), paste(n_trn, "tRNAs"))[which(c(gene.no, rrn.no, trn.no))]
text <- c(paste(ns, collapse = "; "), paste(prettyNum(l, big.mark = ","), "bp", " "),
paste("GC:",round(gc_total, 2)*100,"%"))
j <- 1
if (cu.bias){
cex = 0.7 * text.size
} else {
cex = 0.8 * text.size
}
for (i in which(c((length(ns) != 0), genome.length, total.gc))) {
graphics::text(0,0.05 - 0.05 * j, text[i],
font=4, cex = cex, col = info.color)
j <- j + 1
}
# Highlight IR
# Add legend
if (legend){
legend(x = -1.27, y = -0.6, legend = color_table$label, cex = 0.6 * text.size,
fill = color_table$col, bg = rgb(255, 255, 255, alpha = 0, maxColorValue = 255),
box.col = rgb(255, 255, 255, alpha = 0, maxColorValue = 255))
}
if (save){
grDevices::dev.off()
}
}
shuyuzheng/Chloroplot documentation built on Dec. 23, 2021, 1:28 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions PlotMitGenome PlotPlastidGenome PlotTab
Documented in PlotMitGenome PlotPlastidGenome PlotTab
#' Generate table for PlotGenome function
#'
#' @param gbfile A character. It is the GI or GenBank accession for the
#' interested specie's chloroplast genome files.
#' @param local.file A logical value. If it is \code{TRUE}, the value passed to
#' \code{gbfile} should be a local GB file.
#' @param gc.window A integer. It indicates the window sized for plot
#' GC count.
#'
#' @return A list. It contains 3 tables:
#' \itemize{
#' \item ir_table: a data frame contains information for IR region.
#' \item gene_table: a data frame contains information of genes.
#' \item gc_count: a data frame contains GC count information.
#' }
#'
#' @export
#' @importFrom magrittr %>%
#' @import dplyr
PlotTab <- function(gbfile, local.file = FALSE, gc.window = 100){
my_env <- new.env(parent = emptyenv())
my_env$gb <- NULL
my_env$def <- NULL
my_env$sp_name <- NULL
my_env$genome <- NULL
my_env$l <- NULL
my_env$gene_table <- NULL
my_env$ir <- NULL
if (local.file){
tryCatch({
my_env$gb <- genbankr::readGenBank(gbfile)
my_env$genome <- genbankr::getSeq(my_env$gb)[[1]]
#genome <- rdnFixer(genome)
my_env$l<- Biostrings::nchar(my_env$genome)
my_env$def <- my_env$gb@definition
my_env$sp_name <- sp.name(my_env$def)
my_env$gene_table <- geneTableRead(my_env$gb, my_env$genome)
}, error = function(e){
my_env$gb <- genbankr::parseGenBank(gbfile)
my_env$def <- my_env$gb$DEFINITION
my_env$sp_name <- sp.name(my_env$def)
my_env$genome <- my_env$gb$ORIGIN[[1]]
my_env$l <- Biostrings::nchar(my_env$genome)
my_env$gene_table <- geneTableParsed(my_env$gb, my_env$genome)
})
} else {
tryCatch({
my_env$gb <- genbankr::readGenBank(text = fetch.gb(gbfile))
my_env$genome <- genbankr::getSeq(my_env$gb)
my_env$l <- Biostrings::nchar(my_env$genome)
my_env$def <- my_env$gb@definition
my_env$sp_name <- sp.name(my_env$df)
my_env$gene_table <- geneTableRead(my_env$gb, my_env$genome)
}, error = function(e){
my_env$gb <- genbankr::parseGenBank(text = fetch.gb(gbfile))
my_env$def <- my_env$gb$DEFINITION
my_env$sp_name <- sp.name(my_env$def)
my_env$genome <- my_env$gb$ORIGIN[[1]]
my_env$l <- Biostrings::nchar(my_env$genome)
my_env$gene_table <- geneTableParsed(my_env$gb, my_env$genome)
})
}
# 1. Organelle ------------------------------------------------------------
plastid <- TRUE
if (grepl("mitochondrion", my_env$def)){
plastid <- FALSE
}
# 2. IR LSC SSC -----------------------------------------------------------
# when seed.size is too short it will raise some error in irDetect function
# If it is too long, it cannot correctly detect the IR resions has mutation
# closing the ends.
if (plastid){
tryCatch({
my_env$ir <- irDetect(my_env$genome, seed.size = 100)
}, error = function(e){
my_env$ir <- irDetect(my_env$genome, seed.size = 1000)
}, warning = function(w) {
my_env$ir <- irDetect(my_env$genome, seed.size = 1000)
})
gene_class <- c("psa","psb","pet","atp","ndh","rbc","rpo","rps","rpl",
"clp|mat|inf","ycf","trn","rrn")
} else {
my_env$ir <- list(ir_table = NULL, indel_table = NULL)
gene_class <- c("nad|nd","sdh","cob","cox|cytb","atp","ccmF","mtt","rps","rpl",
"mat","orf","trn","rrn")
}
my_env$gene_table$class <- rep("OTHER", nrow(my_env$gene_table))
for (i in gene_class){
my_env$gene_table$class[which(grepl(as.character(i),
my_env$gene_table$gene, perl = TRUE,
ignore.case = TRUE))] <- i
}
# 3. GC count -------------------------------------------------------------
if (my_env$l > 500000){
gc.window <- 200
} else if (my_env$l < 100000){
gc.window <- 50
}
gc_count_list <- gc_count(my_env$genome, view.width = gc.window)
gc_count <- gc_count_list[[1]]
gc_total <- gc_count_list[[2]]
gc_count$chr <- rep("chr1", nrow(gc_count))
gc_count <- select(gc_count, chr, position, gc_count)
tables <- list(ir_table = my_env$ir$ir_table, indel_table = my_env$ir$indel_table,
gc_count = gc_count, gc_total = gc_total, plastid = plastid,
sp_name = my_env$sp_name, genome_len = my_env$l, gene_table = my_env$gene_table,
genome = my_env$genome)
return(tables)
}
#' Generate Plasit Genome Plot
#'
#' @param plot.tables a list contains information of IR region, gene, and gc
#' count of the genome. It can be generated by function \code{\link{PlotTab}}.
#' @param save A logical value. If it is \code{TRUE}, the plot will be saved in
#' work directory.
#' @param file.type A charactor. It indicates the format of the file in which
#' the plot will be saved. Options are: "pdf", "png", "jpeg","bmp", "tiff".
#' By default, it is set as "pdf".
#' @param file.name A charactor. It indicates the name of the file in which
#' the plot will be saved. By default, it is set as the specie's name.
#' @param shadow A logical value. If it is \code{TRUE}, the section of IR
#' regions will be highlighted by shadows.
#' @param shadow.color An R color object. It indicates the color for the shadow
#' casted from IR sectors.
#' @param legend A logical value. If it is \code{TRUE}, the legend for gene
#' colors will be shown.
#' @param show.indel A logical value. If it is \code{TRUE}, the SNP, insertion
#' and deletion areas in IR regions will be highlighted on the plot.
#' @param ssc.converse A logical value. If it is \code{TRUE}, the SSC region
#' will be converted to its reverse complementary version
#' @param lsc.converse A logical value. If it is \code{TRUE}, the LSC region
#' will be converted to its reverse complementary version
#' @param ira.converse A logical value. If it is \code{TRUE}, the IRA region
#' will be converted to its reverse complementary version
#' @param irb.converse A logical value. If it is \code{TRUE}, the IRB region
#' will be converted to its reverse complementary version
#' @param background An R color object. It indicates the color for the background of
#' entire plot area.
#' @param height A vector of numeric value. The elements of it indicat the
#' height of gene layer, GC count layer and IR region layer, respectively.
#' Default setting is "0.1, 0.2, 0.07". The taltal circle plot region always
#' has a radius of 1, so a height of 0.1 means 10\% of the circle radius.
#' @param gc.color An R color object. It indicates the color for the lines in gc count
#' plot.
#' @param gc.background An R color object. It indicates the color for the background
#' of gc count plot area.
#' @param info.background An R color object. It indicates the color for the background
#' of central area where species' information was shown.
#' @param ir.color An R color object. It indicates the color for the background
#' of IR sectors.
#' @param ssc.color An R color object. It indicates the color for the background
#' of SSC sectors.
#' @param lsc.color An R color object. It indicates the color for the background
#' of LSC sectors.
#' @param ir.gc A logical value. If it is \code{TRUE}, the gc content for
#' IR, LSC, and SSC regions will be plot with deep colors.
#' @param gc.per.gene A logical value. If it is \code{TRUE}, the gc content for
#' each gene will be plot on the gene layer with deep colors.
#' @param pseudo A logical value. If it is \code{TRUE}, the pseudo genes (if
#' there is some in the species) will be marked with a "*" at the end of the
#' labels.
#' @param cu.bias A logical value. If it is \code{TRUE}, the condon usage bias
#' metric for each gene will be shown in the labels. The metric “Measure
#' Independent of Length and Composition (MILC)” was used for evaluate the bias.
#' You can find more details in the reference paper.
#' @param text.size A numeric value. It indicates the size of all texts in the
#' plot. For exmple, \code{text.size = 1.5} means all the texts in the plot will
#' be enlarged as 1.5 times of their original size.
#' @param gc.per.gene A logical value. If it is \code{TRUE}, the GC content of
#' each gene will be show by a darker part in gene rectangles.
#' @param genome.length A logical value. If it is \code{TRUE}, the length of
#' genome will be shown in the center of the plot.
#' @param total.gc A logical value. If it is \code{TRUE}, the GC content of
#' whole genome will be shown in the center of the plot.
#' @param gene.no A logical value. If it is \code{TRUE}, the number of genes in
#' the genome will be shown in the center of the plot.
#' @param rrn.no A logical value. If it is \code{TRUE}, the number of rRNAs in
#' the genome will be shown in the center of the plot.
#' @param trn.no A logical value. If it is \code{TRUE}, the number of tRNAs in
#' the genome will be shown in the center of the plot.
#' @param organelle A logical value. If it is \code{TRUE}, the organelle type of
#' the genome will be shown in the center of the plot.
#' @param psa.color An R color object. It indicates the color for genes of
#' photosystem I
#' @param psb.color An R color object. It indicates the color for genes of
#' photosystem II
#' @param pet.color An R color object. It indicates the color for genes of
#' cytochrome b/f complex
#' @param atp.color An R color object. It indicates the color for genes of
#' ATP synthesis
#' @param ndh.color An R color object. It indicates the color for genes of
#' NADH dehydrogenase
#' @param rbc.color An R color object. It indicates the color for genes of
#' RubisCO larg subunit
#' @param rpo.color An R color object. It indicates the color for genes of
#' RNA polymerase
#' @param rsp.color An R color object. It indicates the color for genes of
#' small ribosomal protein
#' @param rpl.color An R color object. It indicates the color for genes of
#' large ribosomal protein
#' @param clp_mat_inf.color An R color object. It indicates the color for genes
#' of clpP, matK, or infA
#' @param ycf.color An R color object. It indicates the color for genes of
#' hypothetical reading frame
#' @param trn.color An R color object. It indicates the color for genes of
#' transfer RNA
#' @param rrn.color An R color object. It indicates the color for genes of
#' ribosomal RNA
#' @param other_gene.color An R color object. It indicates the color for other
#' genes
#' @param show.gene A vector of characters. It indicates which classes of genes
#' will be shown on the plot. A valiable values are "psa","psb","pet","atp",
#' "ndh","rbc", "rpo","rps","rpl", "clp|mat|inf","ycf", "trn","rrn", "OTHER"
#' @param gene_axis_ir.color An R color object. It indicates the color for other
#' genes.
#' @param gene_axis_lsc.color An R color object. It indicates the color for other
#' genes.
#' @param gene_axis_ssc.color An R color object. It indicates the color for other
#' genes.
#' @param customize.ring1 A data frame. It must contain 2 columns:
#' \itemize{
#' \item \strong{position}: 1-base genomic coordinate for the features.
#' \item \strong{value}: the values for the features.
#' }
#'
#' @param customize.ring1.type A charactor. It indicate the plot type in
#' customize.ring1. Avaliable values are "line", "line + filling", "line + dot",
#' "line + dot + filling", "step line", "step line + filling", "vertical line"
#' @param customize.ring1.color An R color object. It indicates the color for
#' the plots in customized ring 1.
#' @param customize.ring2 A data frame. It must contain 2 columns:
#' \itemize{
#' \item \strong{position}: 1-base genomic coordinate for the features.
#' \item \strong{value}: the values for the features.
#' }
#' @param customize.ring2.type A charactor. It indicate the plot type in
#' customize.ring2. Avaliable values are "line", "line + filling", "line + dot",
#' "line + dot + filling", "step line", "step line + filling", "vertical line"
#' @param customize.ring2.color An R color object. It indicates the color for
#' the plots in customized ring 2.
#' @param customize.ring3 A data frame. It must contain 2 columns:
#' \itemize{
#' \item \strong{start}: 1-base genomic coordinate for the start point of the
#' features.
#' \item \strong{end}: 1-base genomic coordinate for the end point of the
#' features.
#' \item \strong{value}: the values for the features.
#' }
#' @param customize.ring3.color An R color object. It indicates the color for
#' the plots in customized ring 3.
#' @return A plot for chloroplast genome.
#'
#' @references Supek, Fran, and Kristian Vlahovicek. “Comparison of codon usage
#' measures and their applicability in prediction of microbial gene
#' expressivity.” BMC bioinformatics vol. 6 182. 19 Jul.
#' 2005, doi:10.1186/1471-2105-6-182
#'
#' @importFrom magrittr %>%
#' @import dplyr
#' @export
PlotPlastidGenome <- function(plot.tables, save = TRUE, file.type = "pdf",
text.size = 1, height = c(0.1, 0.2, 0.07),
file.name = NULL, shadow = TRUE, ir.gc = TRUE,
gc.per.gene = TRUE, pseudo = TRUE, legend = TRUE,
ssc.converse = FALSE, lsc.converse = FALSE,
ira.converse = FALSE, irb.converse = FALSE,
genome.length = TRUE, total.gc = TRUE, show.indel = TRUE,
gene.no = TRUE, rrn.no = TRUE,trn.no = TRUE,
organelle_type = TRUE,
background = "grey90",gc.color = "grey30",
gc.background = "grey70", info.background = "black",
ir.color = "#2F3941", ssc.color = "#82B6E2",
lsc.color = "#299E96", shadow.color = "#0000FF20",
psa.color = "#2A6332", psb.color = "#4C8805",
pet.color = "#7F992C", atp.color = "#9FBB3D",
ndh.color = "#FEEE50", rbc.color = "#4D9E3F",
rpo.color = "#AE2D29", rsp.color = "#D6AD7C",
rpl.color = "#9C7A4B", clp_mat_inf.color = "#D9662D",
ycf.color = "#71B8A9", trn.color = "#172C7F",
rrn.color = "#D1382A", other_gene.color = "#7D7D7D",
show.genes = c("psa","psb","pet","atp","ndh","rbc",
"rpo","rps","rpl", "clp|mat|inf","ycf",
"trn","rrn", "OTHER"),
gene_axis_ir.color = ir.color,
gene_axis_ssc.color = ssc.color,
gene_axis_lsc.color = lsc.color,
cu.bias = TRUE, customize.ring1 = NULL,
customize.ring1.type = "line",
customize.ring1.color = gc.color,
customize.ring2 = NULL,
customize.ring2.type = "line",
customize.ring2.color = gc.color,
customize.ring3 = NULL,
customize.ring3.color = gc.color){
# Unpack plot.table
genome <- plot.tables$genome
ir_table <- plot.tables$ir_table
gc_count <- plot.tables$gc_count
gc_total <- plot.tables$gc_total
sp_name <- plot.tables$sp_name
gene_table <- plot.tables$gene_table
genome <- plot.tables$genome
indel_table <- plot.tables$indel_table
l <- Biostrings::nchar(genome)
# 1. Modify gene table
# ssc covert
if (ssc.converse){
if (nrow(ir_table) == 1){
warning("Didn't get IR region from thid species. ",
"It's impossible to convert SSC region.")
} else{
tmp <- convert_region(ir_table = ir_table, l = l, region = "SSC",
genome = genome, gene_table = gene_table,
customize.ring1 = customize.ring1,
customize.ring2 = customize.ring2,
customize.ring3 = customize.ring3)
gene_table <- tmp$gene_table
gc_count <- tmp$gc_count
gc_total <- tmp$gc_total
customize.ring1 <- tmp$customize.ring1
customize.ring2 <- tmp$customize.ring2
customize.ring3 <- tmp$customize.ring3
}
}
# LSC covert
if (lsc.converse){
if (nrow(ir_table) == 1){
warning("Didn't get IR region from thid species. ",
"It's impossible to convert LSC region.")
gene_table <- gene_table
} else{
tmp <- convert_region(ir_table = ir_table, l = l, region = "LSC",
genome = genome, gene_table = gene_table,
customize.ring1 = customize.ring1,
customize.ring2 = customize.ring2,
customize.ring3 = customize.ring3)
gene_table <- tmp$gene_table
gc_count <- tmp$gc_count
gc_total <- tmp$gc_total
customize.ring1 <- tmp$customize.ring1
customize.ring2 <- tmp$customize.ring2
customize.ring3 <- tmp$customize.ring3
}
}
# IRA covert
if (ira.converse){
if (nrow(ir_table) == 1){
warning("Didn't get IR region from thid species. ",
"It's impossible to convert IRA region.")
} else{
tmp <- convert_region(ir_table = ir_table, l = l, region = "IRA",
genome = genome, gene_table = gene_table,
customize.ring1 = customize.ring1,
customize.ring2 = customize.ring2,
customize.ring3 = customize.ring3,
indel_table = indel_table)
gene_table <- tmp$gene_table
gc_count <- tmp$gc_count
gc_total <- tmp$gc_total
customize.ring1 <- tmp$customize.ring1
customize.ring2 <- tmp$customize.ring2
customize.ring3 <- tmp$customize.ring3
indel_table <- tmp$indel_table
}
}
# IRB covert
if (irb.converse){
if (nrow(ir_table) == 1){
warning("Didn't get IR region from thid species. ",
"It's impossible to convert IRB region.")
gene_table <- gene_table
} else{
tmp <- convert_region(ir_table = ir_table, l = l, region = "IRB",
genome = genome, gene_table = gene_table,
customize.ring1 = customize.ring1,
customize.ring2 = customize.ring2,
customize.ring3 = customize.ring3,
indel_table = indel_table)
gene_table <- tmp$gene_table
gc_count <- tmp$gc_count
gc_total <- tmp$gc_total
customize.ring1 <- tmp$customize.ring1
customize.ring2 <- tmp$customize.ring2
customize.ring3 <- tmp$customize.ring3
indel_table <- tmp$indel_table
}
}
# colouring
color_table <- plastidGeneColor(psa.color = psa.color, psb.color = psb.color,
pet.color = pet.color, atp.color = atp.color,
ndh.color = ndh.color, rbc.color = rbc.color,
rpo.color = rpo.color, rsp.color = rsp.color,
rpl.color = rpl.color, clp_mat_inf.color = clp_mat_inf.color,
ycf.color = ycf.color, trn.color = trn.color,
rrn.color = rrn.color, other_gene.color = other_gene.color)
gene_table <- dplyr::left_join(gene_table, color_table, by = c("class" = "acronym"))
gene_table <- gene_table[which(gene_table$class %in% show.genes), ]
gene_table$gene <- gsub('(.{1,20})(\\s)', '\\1\n', gene_table$gene)
color_table <- color_table[which(color_table$acronym %in% unique(gene_table$class)), ]
# gene labels
gene_table$gene_label <- gene_table$gene
if (cu.bias){
gene_table$gene_label <- apply(gene_table[, c("gene_label", "cu_bias")], 1,
function(x){
if (is.na(x[2])){
return(x[1])
} else {
return(paste0(x[1], " (",
signif(as.numeric(x[2]), 2), ")"))
}
})
}
if (pseudo){
gene_table$gene_label <- apply(gene_table[, c("gene_label", "pseudo")], 1,
function(x){
if (x[2]){
return(paste0(x[1], " *"))
} else {
return(x[1])
}
})
}
genome <- genome
gene_table_r <- filter(gene_table, strand == "-") %>%
select(chr, start, end, gene, col, gc_count, gene_label) %>%
arrange(start)
gene_table_f <- filter(gene_table, strand == "+") %>%
select(chr, start, end, gene, col, gc_count, gene_label) %>%
arrange(start)
# 2. Set colors
ir_table$bg_col <- rep(NA, nrow(ir_table))
ir_table$bg_col[ir_table$name == "LSC"] <- lsc.color
ir_table$bg_col[ir_table$name == "SSC"] <- ssc.color
ir_table$bg_col[grepl("IR", ir_table$name)] <- ir.color
ir_table$bg_col[ir_table$name == "Genome"] <- lsc.color
ir_table$axis_col <- rep(NA, nrow(ir_table))
ir_table$axis_col[ir_table$name == "LSC"] <- gene_axis_lsc.color
ir_table$axis_col[ir_table$name == "SSC"] <- gene_axis_ssc.color
ir_table$axis_col[grepl("IR", ir_table$name)] <- gene_axis_ir.color
ir_table$axis_col[ir_table$name == "Genome"] <- gene_axis_lsc.color
# Automatically adjust colors
info.color <- CompColor(info.background)
ir_table$inf_col <- rep(NA, nrow(ir_table))
ir_table$inf_col[ir_table$name=="LSC"] <- CompColor(lsc.color)
ir_table$inf_col[ir_table$name=="SSC"] <- CompColor(ssc.color)
ir_table$inf_col[grepl("IR", ir_table$name)] <- CompColor(ir.color)
ir_table$inf_col[ir_table$name == "Genome"] <- CompColor(lsc.color)
if (!is.null(indel_table)){
indel_table$bg_color <- rep("White", #Transparent("#FFFFFF", 0.9),
nrow(indel_table))
}
ir_table$y <- rep(0.25, nrow(ir_table))
ir_table$y[grepl("IR", ir_table$name)] <- 0.75
# 3. Genome information
n_trn <- sum(grepl("trn", gene_table$gene, perl = TRUE))
n_rrn <- sum(grepl("rrn", gene_table$gene, perl = TRUE))
n_gene <- nrow(gene_table) - n_trn - n_rrn
# 4. Generate plot
# Initialize the plot device
if (save){
if (is.null(file.name)){
file <- sp_name
} else {
file <- file.name
}
if (!file.type %in% c("jpeg", "bmp", "png", "tiff", "pdf")){
warning("Can not save plot in ", file.type, " format. Avaliable formats
are 'jpeg', 'bmp', 'png', 'tiff',and 'pdf'.")
} else if (file.type == "pdf") {
grDevices::pdf(paste(file, file.type, sep="."), width=10.3, height=8.9)
} else {
do.call(file.type, args=list(filename=paste(file, file.type, sep="."),
width = 10.3, height = 8.9, units = "in",
res = 1000))
}
}
# Initialize the layout
ssc_deg<-ir_table$center[which(ir_table$name=="SSC")][1]
if ((length(ssc_deg) == 1) & !is.na(ssc_deg)){
ssc_deg <- 360 - round(((ssc_deg)/l) * 360)
rotate <- 270 - ssc_deg
} else {
rotate <- -90
}
# Track indexes
track_index <- 1:8
names(track_index) <- c("gene out 1", "gene out 2", "gene", "gene inner 1",
"gene inner 2", "arrow", "gc count", "ir region")
if (!is.null(customize.ring1)){
track_index[c("arrow", "gc count", "ir region")] <- track_index[c("arrow", "gc count", "ir region")] + 1
track_index["customize 1"] <- track_index["gene inner 2"] + 1
}
if (!is.null(customize.ring2)){
track_index[c("arrow", "gc count", "ir region")] <- track_index[c("arrow", "gc count", "ir region")] + 1
track_index["customize 2"] <- track_index["arrow"] - 1
}
if (!is.null(customize.ring3)){
track_index[c("arrow", "gc count", "ir region")] <- track_index[c("arrow", "gc count", "ir region")] + 1
track_index["customize 3"] <- track_index["arrow"] - 1
}
# Track heights
n_customize_ring <- sum(c(!is.null(customize.ring1), !is.null(customize.ring2), !is.null(customize.ring3)))
if (n_customize_ring == 1){
height[2] <- height[2]/1.5
} else if (n_customize_ring == 2) {
height[2] <- height[2]/2
height[3] <- height[3] * 0.7
} else if (n_customize_ring == 3) {
height[2] <- height[2]/2
height[3] <- height[3] * 0.7
height[1] <- height[1] *0.7
}
circlize::circos.clear()
circlize::circos.par(start.degree = rotate,
gap.after = 0,
track.margin = c(0, 0), cell.padding = c(0, 0, 0, 0))
circlize::circos.genomicInitialize(data=data.frame(chr="chr1", start=0, end=l,
stringsAsFactors = FALSE),
plotType = NULL)
# 1.2. gene label outside
if (nrow(gene_table_r) == 0 ) {
# Background
circlize::draw.sector(0, 360,
rou1 = 1,
rou2 = 0,
col = background, border = NA)
} else {
circlize::circos.genomicLabels(gene_table_r, labels = gene_table_r$gene_label,
side = "outside",cex = 0.5 * text.size,
connection_height = circlize::convert_height(3, "mm"),
labels_height = max(graphics::strwidth(gene_table_r$gene_label,
cex = 0.5 * text.size,
font = graphics::par("font")))
)
# Background
circlize::draw.sector(0, 360,
rou1 = circlize::get.cell.meta.data("cell.bottom.radius",
track.index = 2),
rou2 = 0,
col = background, border = NA)
}
# 3. Gene rectangles
if (nrow(gene_table_r) != 0 & nrow(gene_table_f) != 0) {
circlize::circos.genomicTrack(gene_table_f, factors = as.factor("chr1"),
ylim = c(-1, 1), bg.border = NA,
track.height = height[1], #circlize::convert_height(10, "mm"),
panel.fun = function(region, value, ...) {
if(gc.per.gene){
# genes lay on forward chain
circlize::circos.genomicRect(
gene_table_f[, c("start", "end")],
value = gene_table_f$gene,
ybottom = - gene_table_f$gc_count,
ytop = 0,
col = gene_table_f$col,
border = NA)
circlize::circos.genomicRect(
gene_table_f[, c("start", "end")],
value = gene_table_f$gene,
ybottom = -1,
ytop = 0,
lwd = 0.5,
col = Transparent(gene_table_f$col, 0.7),
border = NA
)
# genes lay on reverse chain
circlize::circos.genomicRect(
gene_table_r[, c("start", "end")],
value = gene_table_r$gene,
ybottom = 0,
ytop = gene_table_r$gc_count,
col = gene_table_r$col,
border = NA)
circlize::circos.genomicRect(
gene_table_r[, c("start", "end")],
value = gene_table_r$gene,
ybottom = 0,
ytop = 1,
lwd = 0.5,
col = Transparent(gene_table_r$col, 0.7),
border = NA
)
} else {
# genes lay on forward chain
circlize::circos.genomicRect(
gene_table_f[, c("start", "end")],
value = gene_table_f$gene,
ybottom = -1,
ytop = 0,
lwd = 0.5,
col = gene_table_f$col
)
circlize::circos.genomicRect(
gene_table_r[, c("start", "end")],
value = gene_table_r$gene,
ybottom = 0,
ytop = 1,
lwd = 0.5,
col = gene_table_r$col
)
}
# the bars in the middle indicate the IR, SSR, LSR
circlize::circos.rect(xleft = ir_table$start,
ybottom = - 0.05,
xright = ir_table$end,
ytop = 0.05,
col = ir_table$axis_col,
border = NA)
})
} else if (nrow(gene_table_f) != 0){
circlize::circos.genomicTrack(gene_table_f, factors = as.factor("chr1"),
ylim = c(-1, 1), bg.border = NA,
track.height = height[1], #circlize::convert_height(10, "mm"),
panel.fun = function(region, value, ...) {
if (gc.per.gene){
# genes lay on forward chain
circlize::circos.genomicRect(
gene_table_f[, c("start", "end")],
value = gene_table_f$gene,
ybottom = - gene_table_f$gc_count,
ytop = 0,
col = gene_table_f$col,
border = NA)
circlize::circos.genomicRect(
gene_table_f[, c("start", "end")],
value = gene_table_f$gene,
ybottom = -1,
ytop = 0,
lwd = 0.5,
col = Transparent(gene_table_f$col, 0.7),
border = NA
)
} else {
# genes lay on forward chain
circlize::circos.genomicRect(
gene_table_f[, c("start", "end")],
value = gene_table_f$gene,
ybottom = -1,
ytop = 0,
lwd = 0.5,
col = gene_table_f$col
)
}
# the bars in the middle indicate the IR, SSR, LSR
circlize::circos.rect(xleft = ir_table$start,
ybottom = - 0.05,
xright = ir_table$end,
ytop = 0.05,
col = ir_table$axis_col,
border = NA)
})
track_index[3:length(track_index)] <- track_index[3:length(track_index)] - 2
track_index <- track_index[-which(names(track_index) %in% c("gene out 1", "gene out 2"))]
} else if (nrow(gene_table_r) != 0) {
circlize::circos.genomicTrack(gene_table_r, factors = as.factor("chr1"),
ylim = c(-1, 1), bg.border = NA,
track.height = height[1], #circlize::convert_height(10, "mm"),
panel.fun = function(region, value, ...) {
if(gc.per.gene){
# genes lay on reverse chain
circlize::circos.genomicRect(
gene_table_r[, c("start", "end")],
value = gene_table_r$gene,
ybottom = 0,
ytop = gene_table_r$gc_count,
col = gene_table_r$col,
border = NA)
circlize::circos.genomicRect(
gene_table_r[, c("start", "end")],
value = gene_table_r$gene,
ybottom = 0,
ytop = 1,
lwd = 0.5,
col = Transparent(gene_table_r$col, 0.7),
border = NA
)
} else {
circlize::circos.genomicRect(
gene_table_r[, c("start", "end")],
value = gene_table_r$gene,
ybottom = 0,
ytop = 1,
lwd = 0.5,
col = gene_table_r$col
)
}
# the bars in the middle indicate the IR, SSR, LSR
circlize::circos.rect(xleft = ir_table$start,
ybottom = - 0.05,
xright = ir_table$end,
ytop = 0.05,
col = ir_table$axis_col,
border = NA)
})
track_index[(which(names(track_index) == "gene inner 2") + 1):length(track_index)] <-
track_index[(which(names(track_index) == "gene inner 2") + 1):length(track_index)] - 2
track_index <- track_index[-which(names(track_index) %in% c("gene inner 1", "gene inner 2"))]
} else {
track_index <- track_index[6:length(track_index)] -5
}
# 4.5. gene labels inside
if (nrow(gene_table_f) != 0){
circlize::circos.genomicLabels(gene_table_f, labels = gene_table_f$gene_label,
side = "inside", cex = 0.5 * text.size,
connection_height = circlize::convert_height(3, "mm"),
labels_height = max(graphics::strwidth(gene_table_f$gene_label,
cex = 0.5 * text.size,
font = graphics::par("font"))))
}
# customized ring 1
if (!is.null(customize.ring1)) {
customize.ring1 <- customize.ring1[order(customize.ring1$position), ]
customize.ring1$chr <- rep("chr1", nrow(customize.ring1))
ylim <- c(min(customize.ring1$value), max(customize.ring1$value))
bar_color <- ifelse(customize.ring1$value > 0, "orangered3", "steelblue")
if (min(customize.ring1$value) > 0) {
ylim[1] <- 0
bar_color <- rep(customize.ring1.color, nrow(customize.ring1))
}
if (max(customize.ring1$value) < 0){
ylim[2] <- 0
bar_color <- rep(customize.ring1.color, nrow(customize.ring1))
}
style <- switch(customize.ring1.type,
"line" = list(type = "l", arear = FALSE, col = customize.ring1.color, lwd = 0.5),
"line + filling" = list(type = "l", arear = TRUE, col = customize.ring1.color, lwd = 0.01),
"line + dot" = list(type = "o", arear = FALSE, col = bar_color, lwd = 0.5),
"line + dot + filling" = list(type = "o", arear = TRUE, col = bar_color, lwd = 0.01),
"step line" = list(type = "s", arear = FALSE, col = customize.ring1.color, lwd = 0.5),
"step line + filling" = list(type = "s", arear = TRUE, col = customize.ring1.color, lwd = 0.01),
"vertical line" = list(type = "h", arear = FALSE, col = bar_color, lwd = 0.5))
circlize::circos.track(factors =as.factor(customize.ring1$chr),
x = customize.ring1$position,
y = customize.ring1$value, ylim = ylim,
track.height = height[2]/2,
track.margin = c(0, circlize::convert_height(1, "mm")),
bg.border = NA, bg.col = gc.background,
panel.fun = function(x, y) {
circlize::circos.lines(x, y, type = style$type,
area = style$area,
baseline = 0,
pch = 1,
cex = 0.2 * text.size,
col = style$col,
lwd = style$lwd)
circlize::circos.segments(x0=0, x1=l, y0=0,
y1=0, lwd=0.5,
lty="16",
col=darken(customize.ring1.color, 0.7))
})
}
# customized ring 2
if (!is.null(customize.ring2)) {
customize.ring2 <- customize.ring2[order(customize.ring2$position), ]
customize.ring2$chr <- rep("chr1", nrow(customize.ring2))
ylim <- c(min(customize.ring2$value), max(customize.ring2$value))
bar_color <- ifelse(customize.ring2$value > 0, "orangered3", "steelblue")
if (min(customize.ring2$value) > 0) {
ylim[1] <- 0
bar_color <- rep(customize.ring2.color, nrow(customize.ring2))
}
if (max(customize.ring2$value) < 0){
ylim[2] <- 0
bar_color <- rep(customize.ring2.color, nrow(customize.ring2))
}
style <- switch(customize.ring2.type,
"line" = list(type = "l", arear = FALSE, col = customize.ring2.color, lwd = 0.5),
"line + filling" = list(type = "l", arear = TRUE, col = customize.ring2.color, lwd = 0.01),
"line + dot" = list(type = "o", arear = FALSE, col = bar_color, lwd = 0.5),
"line + dot + filling" = list(type = "o", arear = TRUE, col = bar_color, lwd = 0.01),
"step line" = list(type = "s", arear = FALSE, col = customize.ring2.color, lwd = 0.5),
"step line + filling" = list(type = "s", arear = TRUE, col = customize.ring2.color, lwd = 0.01),
"vertical line" = list(type = "h", arear = FALSE, col = bar_color, lwd = 0.5))
circlize::circos.track(factors =as.factor(customize.ring2$chr),
x = customize.ring2$position,
y = customize.ring2$value, ylim = ylim,
track.height = height[2]/2,
track.margin = c(0, circlize::convert_height(1, "mm")),
bg.border = NA, bg.col = gc.background,
panel.fun = function(x, y) {
circlize::circos.lines(x, y, type = style$type,
area = style$area,
baseline = 0,
pch = 1,
cex = 0.2 * text.size,
col = style$col,
lwd = style$lwd)
circlize::circos.segments(x0=0, x1=l, y0=0,
y1=0, lwd=0.5,
lty="16",
col=darken(customize.ring2.color, 0.7))
})
}
# customize.ring 3
if (!is.null(customize.ring3)) {
customize.ring3$chr <- rep("chr1", nrow(customize.ring3))
ylim <- c(min(customize.ring3$value), max(customize.ring3$value))
if (min(customize.ring3$value) > 0) {
ylim[1] <- 0
}
if (max(customize.ring3$value) < 0){
ylim[2] <- 0
}
circlize::circos.track(factors =as.factor(customize.ring3$chr),
track.margin = c(0, circlize::convert_height(1, "mm")),
track.height = height[2]/2,
ylim = ylim,
bg.border = NA, bg.col = gc.background,
panel.fun = function(x, y) {
if (min(customize.ring3$value) < 0 &
max(customize.ring3$value) > 0) {
df_large <- customize.ring3[which(customize.ring3$value > 0), ]
df_less <- customize.ring3[which(customize.ring3$value < 0), ]
circlize::circos.rect(xleft = df_large$start,
xright = df_large$end,
ybottom = 0,
ytop = df_large$value,
col = "orangered3",
border = NA)
circlize::circos.rect(xleft = df_less$start,
xright = df_less$end,
ybottom = df_less$value,
ytop = 0,
col = "steelblue",
border = NA)
} else if (min(customize.ring3$value) < 0) {
circlize::circos.rect(xleft = customize.ring3$start,
xright = customize.ring3$end,
ybottom = customize.ring3$value,
ytop = 0,
col = customize.ring3.color,
border = NA)
} else {
circlize::circos.rect(xleft = customize.ring3$start,
xright = customize.ring3$end,
ybottom = 0,
ytop = customize.ring3$value,
col = customize.ring3.color,
border = NA)
}
})
}
# 6. Arrow outside gnome axis
circlize::circos.track(ylim = c(0, 1), track.height = 0.06, bg.border = NA,
track.margin = c(0, circlize::convert_height(1, "mm")),
panel.fun = function(x, y) {
circlize::circos.arrow(x1 = l - l %/% 40,
x2 = l,
y = 0.85,
arrow.position="start",
col = MildColor(background),
border = NA,
arrow.head.length = circlize::convert_x(3, "mm"),
width = 0.15)
# circlize::circos.genomicAxis(h = "bottom")
})
# 7. GC content
circlize::circos.track(factors =as.factor(gc_count$chr),
x=gc_count$position,
y = 1- gc_count$gc_count,
ylim = c(0, 1),
track.height = height[2],
bg.border = NA, bg.col = gc.background,
panel.fun = function(x, y) {
circlize::circos.lines(x, y, type = "s", area = TRUE,
baseline = "top",
col = gc.color, lwd = 0.00001)
circlize::circos.arrow(x1= 0,
x2= l %/% 35,
y= 0.9,
col=MildColor(gc.color),
border=NA,
arrow.head.length=circlize::convert_x(3, "mm"),
width=0.15 * 0.05 / height[2])
circlize::circos.segments(x0=0, x1=l, y0=0.25,
y1=0.25, lwd=0.5,
lty="16", col=MildColor(gc.background))
circlize::circos.segments(x0=0, x1=l, y0=0.5,
y1=0.5, lwd=0.6,
lty="11", col=MildColor(gc.background))
circlize::circos.segments(x0=0, x1=l, y0=0.75,
y1=0.75, lwd=0.5,
lty="16", col=MildColor(gc.color))
circlize::circos.genomicAxis(h = "top", labels.cex = 0.4 * text.size)
})
#circos.yaxis(at=c(0.25,0.5,0.75),labels.cex=0.25 * text.size,lwd=0,tick.length=0,
# labels.col="grey40",col="#FFFFFF")
# 8. inner ring for IR, SSR, LSR
circlize::circos.genomicTrack(ir_table, bg.border = NA,
ylim = c(0, 1),track.height = height[3],
panel.fun = function(region, value, ...) {
if (ir.gc){
circlize::circos.rect(xleft = ir_table$start,
ybottom = 0,#circlize::convert_y(0, "mm"),
xright = ir_table$end,
ytop = 1,#circlize::convert_y(3, "mm"),
col = Transparent(ir_table$bg_col, 0.7),
border = NA)
circlize::circos.rect(xleft = ir_table$start,
ybottom =0,#circlize::convert_y(0, "mm"),
xright = ir_table$end,
ytop = ir_table$gc_count,#circlize::convert_y(3, "mm"),
col = ir_table$bg_col,
border = NA)
} else {
circlize::circos.rect(xleft = ir_table$start,
ybottom = 0,#circlize::convert_y(0, "mm"),
xright = ir_table$end,
ytop = 1,#circlize::convert_y(3, "mm"),
col = ir_table$bg_col,
border = NA)
}
if (!is.null(indel_table) & show.indel){
circlize::circos.rect(xleft = indel_table$position - 10,
ybottom = 0,#circlize::convert_y(0, "mm"),
xright = indel_table$position + 10,
ytop = 1,#circlize::convert_y(3, "mm"),
col = indel_table$bg_col,
border = NA)
circlize::circos.points(x = indel_table$position,
y = 0.05,
col = indel_table$col,
pch = 19,
cex = 0.5 * text.size)
# circlize::circos.text(x = indel_table$position,
# y = 0.5,
# col = indel_table$col,
# labels = indel_table$string,
# cex = 0.2 * text.size)
}
circlize::circos.text(x = ir_table$center,
y = ir_table$y, cex = 0.7 * text.size,
labels = ir_table$text,
facing = "bending.inside",
col = ir_table$inf_col,
niceFacing = TRUE)
})
# Specie's information inthe central
circlize::draw.sector(0, 360,
rou1 = circlize::get.cell.meta.data("cell.bottom.radius",
track.index = track_index['ir region']),
rou2 = 0,
col = info.background, border = NA)
if (organelle_type){
graphics::text(0,0.10, sp_name, font=4, cex = 0.8 * text.size, col = info.color)
graphics::text(0,0.05, "Plastid Genome", font=4, cex = 0.8 * text.size, col = info.color)
} else{
graphics::text(0,0.05, sp_name, font=4, cex = 0.8 * text.size, col = info.color)
}
ns <- c(paste(n_gene, "genes"), paste(n_rrn, "rRNAs"), paste(n_trn, "tRNAs"))[which(c(gene.no, rrn.no, trn.no))]
text <- c(paste(ns, collapse = "; "), paste(prettyNum(l, big.mark = ","), "bp", " "),
paste("GC:",round(gc_total, 2)*100,"%"))
j <- 1
if (cu.bias){
cex = 0.7 * text.size
} else {
cex = 0.8 * text.size
}
for (i in which(c((length(ns) != 0), genome.length, total.gc))) {
graphics::text(0,0.05 - 0.05 * j, text[i],
font=4, cex = cex, col = info.color)
j <- j + 1
}
# Highlight IR
if (shadow & (nrow(ir_table) >= 4)){
pos_s=circlize::circlize(ir_table$start[which(grepl("IR",
ir_table$name))],
1, sector.index = "chr1", track.index = 1)
pos_e=circlize::circlize(ir_table$end[which(grepl("IR",
ir_table$name))],
1, sector.index = "chr1", track.index = 1)
for (i in 1:nrow(pos_s)) {
if ('gene' %in% names(track_index)){
circlize::draw.sector(pos_s[i, "theta"], pos_e[i, "theta"],
#start.degree = 0, end.degree = 360 - 10,
rou1 = circlize::get.cell.meta.data("cell.top.radius",
track.index = track_index['gene']),
rou2 = circlize::get.cell.meta.data("cell.bottom.radius",
track.index = track_index['gc count']),
col = shadow.color, border = NA)
} else {
circlize::draw.sector(pos_s[i, "theta"], pos_e[i, "theta"],
#start.degree = 0, end.degree = 360 - 10,
rou1 = circlize::get.cell.meta.data("cell.top.radius",
track.index = track_index['arrow']),
rou2 = circlize::get.cell.meta.data("cell.bottom.radius",
track.index = track_index['gc count']),
col = shadow.color, border = NA)
}
}
if (!is.null(indel_table) & show.indel) {
pos_s=circlize::circlize(indel_table$position - 10,
1, sector.index = "chr1", track.index = 1)
pos_e=circlize::circlize(indel_table$position + 10,
1, sector.index = "chr1", track.index = 1)
for (i in 1:nrow(pos_s)) {
if ('gene' %in% names(track_index)){
circlize::draw.sector(pos_s[i, "theta"], pos_e[i, "theta"],
#start.degree = 0, end.degree = 360 - 10,
rou1 = circlize::get.cell.meta.data("cell.top.radius",
track.index = track_index['gene']),
rou2 = circlize::get.cell.meta.data("cell.bottom.radius",
track.index = track_index['gc count']),
col = Transparent("#FFFFFF", 0.5),
border = NA)
} else {
circlize::draw.sector(pos_s[i, "theta"], pos_e[i, "theta"],
#start.degree = 0, end.degree = 360 - 10,
rou1 = circlize::get.cell.meta.data("cell.top.radius",
track.index = track_index['arrow']),
rou2 = circlize::get.cell.meta.data("cell.bottom.radius",
track.index = track_index['gc count']),
col = Transparent("#FFFFFF", 0.9),
border = NA)
}
}
}
}
# Add legend
if (legend){
legend(x = -1.27, y = -0.6, legend = color_table$label, cex = 0.6 * text.size,
fill = color_table$col, box.col = "white")
}
if (save){
grDevices::dev.off()
}
}
#' Generate Mitochondrion Genome Plot
#'
#' @param plot.tables a list contains information of IR region, gene, and gc
#' count of the genome. It can be generated by function \code{\link{PlotTab}}.
#' @param save A logical value. If it is \code{TRUE}, the plot will be saved in
#' work directory.
#' @param file.type A charactor. It indicates the format of the file in which
#' the plot will be saved. Options are: "pdf", "png", "jpeg","bmp", "tiff".
#' By default, it is set as "pdf".
#' @param file.name A charactor. It indicates the name of the file in which
#' the plot will be saved. By default, it is set as the specie's name.
#' @param legend A logical value. If it is \code{TRUE}, the legend for gene
#' colors will be shown.
#' @param background An R color object. It indicates the color for the background of
#' entire plot area.
#' @param height A vector of numeric value. The elements of it indicat the
#' height of gene layer, GC count layer and IR region layer, respectively.
#' Default setting is "0.1, 0.2, 0.07". The taltal circle plot region always
#' has a radius of 1, so a height of 0.1 means 10\% of the circle radius.
#' @param gc.color An R color object. It indicates the color for the lines in gc count
#' plot.
#' @param gc.background An R color object. It indicates the color for the background
#' of gc count plot area.
#' @param info.background An R color object. It indicates the color for the background
#' of central area where species' information was shown.
#' @param gc.per.gene A logical value. If it is \code{TRUE}, the gc content for
#' each gene will be plot on the gene layer with deep colors.
#' @param pseudo A logical value. If it is \code{TRUE}, the pseudo genes (if
#' there is some in the species) will be marked with a "*" at the end of the
#' labels.
#' @param cu.bias A logical value. If it is \code{TRUE}, the condon usage bias
#' metric for each gene will be shown in the labels. The metric “Measure
#' Independent of Length and Composition (MILC)” was used for evaluate the bias.
#' You can find more details in the reference paper.
#' @param text.size A numeric value. It indicates the size of all texts in the
#' plot. For exmple, \code{text.size = 1.5} means all the texts in the plot will
#' be enlarged as 1.5 times of their original size.
#' @param gc.per.gene A logical value. If it is \code{TRUE}, the GC content of
#' each gene will be show by a darker part in gene rectangles.
#' @param genome.length A logical value. If it is \code{TRUE}, the length of
#' genome will be shown in the center of the plot.
#' @param total.gc A logical value. If it is \code{TRUE}, the GC content of
#' whole genome will be shown in the center of the plot.
#' @param gene.no A logical value. If it is \code{TRUE}, the number of genes in
#' the genome will be shown in the center of the plot.
#' @param rrn.no A logical value. If it is \code{TRUE}, the number of rRNAs in
#' the genome will be shown in the center of the plot.
#' @param trn.no A logical value. If it is \code{TRUE}, the number of tRNAs in
#' the genome will be shown in the center of the plot.
#' @param organelle A logical value. If it is \code{TRUE}, the organelle type of
#' the genome will be shown in the center of the plot.
#' @param nad.color An R color string. It indicates the color for genes of
#' complex I (NADH dehydrogenase)
#' @param sdh.color An R color string. It indicates the color for genes of
#' complex II (succinate dehydrogenase)
#' @param cob.color An R color string. It indicates the color for genes of
#' complex III (ubichinol cytochrome reductase)
#' @param cox.color An R color string. It indicates the color for genes of
#' complex IV (cytochrome c oxidase)
#' @param atp.color An R color string. It indicates the color for genes of
#' ATP synthase
#' @param ccmF.color An R color string. It indicates the color for genes of
#' cytochrome c biogenesi
#' @param mmt.color An R color string. It indicates the color for genes of
#' transport membrane protein
#' @param rps.color An R color string. It indicates the color for genes of
#' ribosomal proteins (SSU)
#' @param rpl.color An R color string. It indicates the color for genes of
#' ribosomal proteins (LSU)
#' @param mat.color An R color string. It indicates the color for genes of
#' maturases
#' @param orf.color An R color string. It indicates the color for ORFs
#' @param trn.color An R color object. It indicates the color for genes of
#' transfer RNA
#' @param rrn.color An R color object. It indicates the color for genes of
#' ribosomal RNA
#' @param other_gene.color An R color object. It indicates the color for other
#' genes
#' @param show.gene A vector of characters. It indicates which classes of genes
#' will be shown on the plot. Avaliable values are "nad|nd","sdh","cob",
#' "cox|cytb","atp", "ccmF","mtt","rps","rpl", "mat","orf","trn","rrn", "OTHER"
#' @param customize.ring1 A data frame. It must contain 2 columns:
#' \itemize{
#' \item \strong{position}: 1-base genomic coordinate for the features.
#' \item \strong{value}: the values for the features.
#' }
#'
#' @param customize.ring1.type A character. It indicate the plot type in
#' customize.ring1. Available values are "line", "line + filling", "line + dot",
#' "line + dot + filling", "step line", "step line + filling", "vertical line"
#' @param customize.ring1.color An R color object. It indicates the color for
#' the plots in customized ring 1.
#' @param customize.ring2 A data frame. It must contain 2 columns:
#' \itemize{
#' \item \strong{position}: 1-base genome coordinate for the features.
#' \item \strong{value}: the values for the features.
#' }
#' @param customize.ring2.type A character. It indicate the plot type in
#' customize.ring2. Available values are "line", "line + filling", "line + dot",
#' "line + dot + filling", "step line", "step line + filling", "vertical line"
#' @param customize.ring2.color An R color object. It indicates the color for
#' the plots in customized ring 2.
#' @param customize.ring3 A data frame. It must contain 2 columns:
#' \itemize{
#' \item \strong{start}: 1-base genome coordinate for the start point of the
#' features.
#' \item \strong{end}: 1-base genome coordinate for the end point of the
#' features.
#' \item \strong{value}: the values for the features.
#' }
#' @param customize.ring3.color An R color object. It indicates the color for
#' the plots in customized ring 3.
#' @return A plot for chloroplast genome.
#'
#' @references Supek, Fran, and Kristian Vlahovicek. “Comparison of codon usage
#' measures and their applicability in prediction of microbial gene
#' expressivity.” BMC bioinformatics vol. 6 182. 19 Jul.
#' 2005, doi:10.1186/1471-2105-6-182
#'
#' @importFrom magrittr %>%
#' @import dplyr
#' @export
PlotMitGenome <- function(plot.tables, save = TRUE, file.type = "pdf",
text.size = 1, height = c(0.1, 0.2, 0.07),
file.name = NULL,
gc.per.gene = TRUE, pseudo = TRUE, legend = TRUE,
genome.length = TRUE, total.gc = TRUE,
gene.no = TRUE, rrn.no = TRUE, trn.no = TRUE,
organelle_type = TRUE,
background = "grey90", gc.color = "grey30",
gc.background = "grey70", info.background = "black",
nad.color = "#2A6332", sdh.color = "#4C8805",
cob.color = "#7F992C", cox.color = "#FEEE50",
atp.color = "#9FBB3D", ccmF.color = "#4D9E3F",
mmt.color = "#AE2D29", rps.color = "#D6AD7C",
rpl.color = "#9C7A4B", mat.color = "#D9662D",
orf.color = "#71B8A9", trn.color = "#172C7F",
rrn.color = "#D1382A", other_gene.color = "#7D7D7D",
show.genes = c("nad|nd","sdh","cob","cox|cytb","atp",
"ccmF","mtt","rps","rpl",
"mat","orf","trn","rrn", "OTHER"),
cu.bias = TRUE, customize.ring1 = NULL,
customize.ring1.color = "grey30",
customize.ring1.type = "line",
customize.ring2.color = "grey30",
customize.ring2 = NULL,
customize.ring2.type = "line",
customize.ring3 = NULL,
customize.ring3.color = "grey30"){
# Unpack plot.table
genome <- plot.tables$genome
gc_count <- plot.tables$gc_count
gc_total <- plot.tables$gc_total
sp_name <- plot.tables$sp_name
gene_table <- plot.tables$gene_table
genome <- plot.tables$genome
l <- Biostrings::nchar(genome)
# 1. Modify gene table
# colouring
color_table <- mitGeneColor(nad.color = nad.color, sdh.color = sdh.color,
cob.color = cob.color, cox.color = cox.color,
atp.color = atp.color, ccmF.color = ccmF.color,
mmt.color = mmt.color, rps.color = rps.color,
rpl.color = mmt.color, mat.color = mat.color,
orf.color = orf.color, trn.color = trn.color,
rrn.color = rrn.color,
other_gene.color = other_gene.color)
gene_table <- dplyr::left_join(gene_table, color_table, by = c("class" = "acronym"))
gene_table <- gene_table[which(gene_table$class %in% show.genes), ]
color_table <- color_table[which(color_table$acronym %in% unique(gene_table$class)), ]
# gene labels
gene_table$gene_label <- gene_table$gene
if (cu.bias){
gene_table$gene_label <- apply(gene_table[, c("gene_label", "cu_bias")], 1,
function(x){
if (is.na(x[2])){
return(x[1])
} else {
return(paste0(x[1], " (",
signif(as.numeric(x[2]), 2), ")"))
}
})
}
if (pseudo){
gene_table$gene_label <- apply(gene_table[, c("gene_label", "pseudo")], 1,
function(x){
if (x[2]){
return(paste0(x[1], " *"))
} else {
return(x[1])
}
})
}
genome <- genome
gene_table_r <- filter(gene_table, strand == "-") %>%
select(chr, start, end, gene, col, gc_count, gene_label) %>%
arrange(start)
gene_table_f <- filter(gene_table, strand == "+") %>%
select(chr, start, end, gene, col, gc_count, gene_label) %>%
arrange(start)
# 2. Set colors
# Automatically adjust colors
info.color <- CompColor(info.background)
# 3. Genome information
n_trn <- sum(grepl("trn", gene_table$gene, perl = TRUE))
n_rrn <- sum(grepl("rrn", gene_table$gene, perl = TRUE))
n_gene <- nrow(gene_table) - n_trn - n_rrn
# 4. Generate plot
# Initialize the plot device
if (save){
if (is.null(file.name)){
file <- sp_name
} else {
file <- file.name
}
if (!file.type %in% c("jpeg", "bmp", "png", "tiff", "pdf")){
warning("Can not save plot in ", file.type, " format. Avaliable formats
are 'jpeg', 'bmp', 'png', 'tiff',and 'pdf'.")
} else if (file.type == "pdf") {
grDevices::pdf(paste(file, file.type, sep="."), width=10.3, height=8.9)
} else {
do.call(file.type, args=list(filename=paste(file, file.type, sep="."),
width = 10.3, height = 8.9, units = "in",
res = 1000))
}
}
# Initialize the layout
# Track indexes
track_index <- 1:7
names(track_index) <- c("gene out 1", "gene out 2", "gene", "gene inner 1",
"gene inner 2", "arrow", "gc count")
if (!is.null(customize.ring1)){
track_index[c("arrow", "gc count")] <- track_index[c("arrow", "gc count")] + 1
track_index["customize 1"] <- track_index["gene inner 2"] + 1
}
if (!is.null(customize.ring2)){
track_index[c("arrow", "gc count")] <- track_index[c("arrow", "gc count")] + 1
track_index["customize 2"] <- track_index["arrow"] - 1
}
if (!is.null(customize.ring3)){
track_index[c("arrow", "gc count")] <- track_index[c("arrow", "gc count")] + 1
track_index["customize 3"] <- track_index["arrow"] - 1
}
# Track heights
n_customize_ring <- sum(c(!is.null(customize.ring1), !is.null(customize.ring2), !is.null(customize.ring3)))
if (n_customize_ring == 1){
height[2] <- height[2]/1.5
} else if (n_customize_ring == 2) {
height[2] <- height[2]/2
height[3] <- height[3] * 0.7
} else if (n_customize_ring == 3) {
height[2] <- height[2]/2
height[3] <- height[3] * 0.7
height[1] <- height[1] *0.7
}
circlize::circos.clear()
circlize::circos.par(start.degree = 0, gap.after = 0,
track.margin = c(0, 0), cell.padding = c(0, 0, 0, 0))
circlize::circos.genomicInitialize(data=data.frame(chr="chr1", start=0, end=l,
stringsAsFactors = FALSE),
plotType = NULL)
# 1.2. gene label outside
if (nrow(gene_table_r) == 0 ) {
# Background
circlize::draw.sector(0, 360,
rou1 = 1,
rou2 = 0,
col = background, border = NA)
} else {
circlize::circos.genomicLabels(gene_table_r, labels = gene_table_r$gene_label,
side = "outside",cex = 0.5 * text.size,
connection_height = circlize::convert_height(3, "mm"),
labels_height = max(graphics::strwidth(gene_table_r$gene_label,
cex = 0.5 * text.size,
font = graphics::par("font")))
)
# Background
circlize::draw.sector(0, 360,
rou1 = circlize::get.cell.meta.data("cell.bottom.radius",
track.index = 2),
rou2 = 0,
col = background, border = NA)
}
# 3. Gene rectangles
if (nrow(gene_table_r) != 0 & nrow(gene_table_f) != 0) {
circlize::circos.genomicTrack(gene_table_f, factors = as.factor("chr1"),
ylim = c(-1, 1), bg.border = NA,
track.height = height[1], #circlize::convert_height(10, "mm"),
panel.fun = function(region, value, ...) {
if(gc.per.gene){
# genes lay on forward chain
circlize::circos.genomicRect(
gene_table_f[, c("start", "end")],
value = gene_table_f$gene,
ybottom = - gene_table_f$gc_count,
ytop = 0,
col = gene_table_f$col,
border = NA)
circlize::circos.genomicRect(
gene_table_f[, c("start", "end")],
value = gene_table_f$gene,
ybottom = -1,
ytop = 0,
lwd = 0.5,
col = Transparent(gene_table_f$col, 0.7),
border = NA
)
# genes lay on reverse chain
circlize::circos.genomicRect(
gene_table_r[, c("start", "end")],
value = gene_table_r$gene,
ybottom = 0,
ytop = gene_table_r$gc_count,
col = gene_table_r$col,
border = NA)
circlize::circos.genomicRect(
gene_table_r[, c("start", "end")],
value = gene_table_r$gene,
ybottom = 0,
ytop = 1,
lwd = 0.5,
col = Transparent(gene_table_r$col, 0.7),
border = NA
)
} else {
# genes lay on forward chain
circlize::circos.genomicRect(
gene_table_f[, c("start", "end")],
value = gene_table_f$gene,
ybottom = -1,
ytop = 0,
lwd = 0.5,
col = gene_table_f$col
)
circlize::circos.genomicRect(
gene_table_r[, c("start", "end")],
value = gene_table_r$gene,
ybottom = 0,
ytop = 1,
lwd = 0.5,
col = gene_table_r$col
)
}
# the x axis
circlize::circos.rect(xleft = 0,
ybottom = - 0.05,
xright = l,
ytop = 0.05,
col = "black",
border = NA)
})
} else if (nrow(gene_table_f) != 0){
circlize::circos.genomicTrack(gene_table_f, factors = as.factor("chr1"),
ylim = c(-1, 1), bg.border = NA,
track.height = height[1], #circlize::convert_height(10, "mm"),
panel.fun = function(region, value, ...) {
if (gc.per.gene){
# genes lay on forward chain
circlize::circos.genomicRect(
gene_table_f[, c("start", "end")],
value = gene_table_f$gene,
ybottom = - gene_table_f$gc_count,
ytop = 0,
col = gene_table_f$col,
border = NA)
circlize::circos.genomicRect(
gene_table_f[, c("start", "end")],
value = gene_table_f$gene,
ybottom = -1,
ytop = 0,
lwd = 0.5,
col = Transparent(gene_table_f$col, 0.7),
border = NA
)
} else {
# genes lay on forward chain
circlize::circos.genomicRect(
gene_table_f[, c("start", "end")],
value = gene_table_f$gene,
ybottom = -1,
ytop = 0,
lwd = 0.5,
col = gene_table_f$col
)
}
# the x axis
circlize::circos.rect(xleft = 0,
ybottom = - 0.05,
xright = l,
ytop = 0.05,
col = "black",
border = NA)
})
track_index[3:length(track_index)] <- track_index[3:length(track_index)] - 2
track_index <- track_index[-which(names(track_index) %in% c("gene out 1", "gene out 2"))]
} else if (nrow(gene_table_r) != 0) {
circlize::circos.genomicTrack(gene_table_r, factors = as.factor("chr1"),
ylim = c(-1, 1), bg.border = NA,
track.height = height[1], #circlize::convert_height(10, "mm"),
panel.fun = function(region, value, ...) {
if(gc.per.gene){
# genes lay on reverse chain
circlize::circos.genomicRect(
gene_table_r[, c("start", "end")],
value = gene_table_r$gene,
ybottom = 0,
ytop = gene_table_r$gc_count,
col = gene_table_r$col,
border = NA)
circlize::circos.genomicRect(
gene_table_r[, c("start", "end")],
value = gene_table_r$gene,
ybottom = 0,
ytop = 1,
lwd = 0.5,
col = Transparent(gene_table_r$col, 0.7),
border = NA
)
} else {
circlize::circos.genomicRect(
gene_table_r[, c("start", "end")],
value = gene_table_r$gene,
ybottom = 0,
ytop = 1,
lwd = 0.5,
col = gene_table_r$col
)
}
# the x axis
circlize::circos.rect(xleft = 0,
ybottom = - 0.05,
xright = l,
ytop = 0.05,
col = "black",
border = NA)
})
track_index[(which(names(track_index) == "gene inner 2") + 1):length(track_index)] <-
track_index[(which(names(track_index) == "gene inner 2") + 1):length(track_index)] - 2
track_index <- track_index[-which(names(track_index) %in% c("gene inner 1", "gene inner 2"))]
} else {
track_index <- track_index[6:length(track_index)] -5
}
# 4.5. gene labels inside
if (nrow(gene_table_f) != 0){
circlize::circos.genomicLabels(gene_table_f, labels = gene_table_f$gene_label,
side = "inside", cex = 0.5 * text.size,
connection_height = circlize::convert_height(3, "mm"),
labels_height = max(graphics::strwidth(gene_table_f$gene_label,
cex = 0.5 * text.size,
font = graphics::par("font"))))
}
# customized ring 1
if (!is.null(customize.ring1)) {
customize.ring1 <- customize.ring1[order(customize.ring1$position), ]
customize.ring1$chr <- rep("chr1", nrow(customize.ring1))
ylim <- c(min(customize.ring1$value), max(customize.ring1$value))
bar_color <- ifelse(customize.ring1$value > 0, "orangered3", "steelblue")
if (min(customize.ring1$value) > 0) {
ylim[1] <- 0
bar_color <- rep(customize.ring1.color, nrow(customize.ring1))
}
if (max(customize.ring1$value) < 0){
ylim[2] <- 0
bar_color <- rep(customize.ring1.color, nrow(customize.ring1))
}
style <- switch(customize.ring1.type,
"line" = list(type = "l", arear = FALSE, col = customize.ring1.color, lwd = 0.5),
"line + filling" = list(type = "l", arear = TRUE, col = customize.ring1.color, lwd = 0.01),
"line + dot" = list(type = "o", arear = FALSE, col = bar_color, lwd = 0.5),
"line + dot + filling" = list(type = "o", arear = TRUE, col = bar_color, lwd = 0.01),
"step line" = list(type = "s", arear = FALSE, col = customize.ring1.color, lwd = 0.5),
"step line + filling" = list(type = "s", arear = TRUE, col = customize.ring1.color, lwd = 0.01),
"vertical line" = list(type = "h", arear = FALSE, col = bar_color, lwd = 0.5))
circlize::circos.track(factors =as.factor(customize.ring1$chr),
x = customize.ring1$position,
y = customize.ring1$value, ylim = ylim,
track.height = height[2]/2,
track.margin = c(0, circlize::convert_height(1, "mm")),
bg.border = NA, bg.col = gc.background,
panel.fun = function(x, y) {
circlize::circos.lines(x, y, type = style$type,
area = style$area,
baseline = 0,
pch = 1,
cex = 0.2 * text.size,
col = style$col,
lwd = style$lwd)
circlize::circos.segments(x0=0, x1=l, y0=0,
y1=0, lwd=0.5,
lty="16",
col=darken(customize.ring1.color, 0.7))
})
}
# customized ring 2
if (!is.null(customize.ring2)) {
customize.ring2 <- customize.ring2[order(customize.ring2$position), ]
customize.ring2$chr <- rep("chr1", nrow(customize.ring2))
ylim <- c(min(customize.ring2$value), max(customize.ring2$value))
bar_color <- ifelse(customize.ring2$value > 0, "orangered3", "steelblue")
if (min(customize.ring2$value) > 0) {
ylim[1] <- 0
bar_color <- rep(customize.ring2.color, nrow(customize.ring2))
}
if (max(customize.ring2$value) < 0){
ylim[2] <- 0
bar_color <- rep(customize.ring2.color, nrow(customize.ring2))
}
style <- switch(customize.ring2.type,
"line" = list(type = "l", arear = FALSE, col = customize.ring2.color, lwd = 0.5),
"line + filling" = list(type = "l", arear = TRUE, col = customize.ring2.color, lwd = 0.01),
"line + dot" = list(type = "o", arear = FALSE, col = bar_color, lwd = 0.5),
"line + dot + filling" = list(type = "o", arear = TRUE, col = bar_color, lwd = 0.01),
"step line" = list(type = "s", arear = FALSE, col = customize.ring2.color, lwd = 0.5),
"step line + filling" = list(type = "s", arear = TRUE, col = customize.ring2.color, lwd = 0.01),
"vertical line" = list(type = "h", arear = FALSE, col = bar_color, lwd = 0.5))
circlize::circos.track(factors =as.factor(customize.ring2$chr),
x = customize.ring2$position,
y = customize.ring2$value, ylim = ylim,
track.height = height[2]/2,
track.margin = c(0, circlize::convert_height(1, "mm")),
bg.border = NA, bg.col = gc.background,
panel.fun = function(x, y) {
circlize::circos.lines(x, y, type = style$type,
area = style$area,
baseline = 0,
pch = 1,
cex = 0.2 * text.size,
col = style$col,
lwd = style$lwd)
circlize::circos.segments(x0=0, x1=l, y0=0,
y1=0, lwd=0.5,
lty="16",
col=darken(customize.ring2.color, 0.7))
})
}
# customize.ring 3
if (!is.null(customize.ring3)) {
customize.ring3$chr <- rep("chr1", nrow(customize.ring3))
ylim <- c(min(customize.ring3$value), max(customize.ring3$value))
if (min(customize.ring3$value) > 0) {
ylim[1] <- 0
}
if (max(customize.ring3$value) < 0){
ylim[2] <- 0
}
circlize::circos.track(factors =as.factor(customize.ring3$chr),
track.margin = c(0, circlize::convert_height(1, "mm")),
track.height = height[2]/2,
ylim = ylim,
bg.border = NA, bg.col = gc.background,
panel.fun = function(x, y) {
if (min(customize.ring3$value) < 0 &
max(customize.ring3$value) > 0) {
df_large <- customize.ring3[which(customize.ring3$value > 0), ]
df_less <- customize.ring3[which(customize.ring3$value < 0), ]
circlize::circos.rect(xleft = df_large$start,
xright = df_large$end,
ybottom = 0,
ytop = df_large$value,
col = "orangered3",
border = NA)
circlize::circos.rect(xleft = df_less$start,
xright = df_less$end,
ybottom = df_less$value,
ytop = 0,
col = "steelblue",
border = NA)
} else if (min(customize.ring3$value) < 0) {
circlize::circos.rect(xleft = customize.ring3$start,
xright = customize.ring3$end,
ybottom = customize.ring3$value,
ytop = 0,
col = customize.ring3.color,
border = NA)
} else {
circlize::circos.rect(xleft = customize.ring3$start,
xright = customize.ring3$end,
ybottom = 0,
ytop = customize.ring3$value,
col = customize.ring3.color,
border = NA)
}
})
}
# 6. Arrow outside gnome axis
circlize::circos.track(ylim = c(0, 1), track.height = 0.06, bg.border = NA,
track.margin = c(0, circlize::convert_height(1, "mm")),
panel.fun = function(x, y) {
circlize::circos.arrow(x1 = l - l %/% 40,
x2 = l,
y = 0.85,
arrow.position="start",
col = MildColor(background),
border = NA,
arrow.head.length = circlize::convert_x(3, "mm"),
width = 0.15)
# circlize::circos.genomicAxis(h = "bottom")
})
# 7. GC content
circlize::circos.track(factors =as.factor(gc_count$chr),
x=gc_count$position,
y = 1- gc_count$gc_count,
ylim = c(0, 1),
track.height = height[2],
bg.border = NA, bg.col = gc.background,
panel.fun = function(x, y) {
circlize::circos.lines(x, y, type = "s", area = TRUE,
baseline = "top",
col = gc.color, lwd = 0.00001)
circlize::circos.arrow(x1= 0,
x2= l %/% 35,
y= 0.9,
col=MildColor(gc.color),
border=NA,
arrow.head.length=circlize::convert_x(3, "mm"),
width=0.15 * 0.05 / height[2])
circlize::circos.segments(x0=0, x1=l, y0=0.25,
y1=0.25, lwd=0.5,
lty="16", col=MildColor(gc.background))
circlize::circos.segments(x0=0, x1=l, y0=0.5,
y1=0.5, lwd=0.6,
lty="11", col=MildColor(gc.background))
circlize::circos.segments(x0=0, x1=l, y0=0.75,
y1=0.75, lwd=0.5,
lty="16", col=MildColor(gc.color))
circlize::circos.genomicAxis(h = "top", labels.cex = 0.4 * text.size)
})
#circos.yaxis(at=c(0.25,0.5,0.75),labels.cex=0.25 * text.size,lwd=0,tick.length=0,
# labels.col="grey40",col="#FFFFFF")
# Specie's information inthe central
circlize::draw.sector(0, 360,
rou1 = circlize::get.cell.meta.data("cell.bottom.radius",
track.index = track_index['gc count']),
rou2 = 0,
col = info.background, border = NA)
if (organelle_type){
graphics::text(0,0.10, sp_name, font=4, cex = 0.8 * text.size, col = info.color)
graphics::text(0,0.05, "Plastid Genome", font=4, cex = 0.8 * text.size, col = info.color)
} else{
graphics::text(0,0.05, sp_name, font=4, cex = 0.8 * text.size, col = info.color)
}
ns <- c(paste(n_gene, "genes"), paste(n_rrn, "rRNAs"), paste(n_trn, "tRNAs"))[which(c(gene.no, rrn.no, trn.no))]
text <- c(paste(ns, collapse = "; "), paste(prettyNum(l, big.mark = ","), "bp", " "),
paste("GC:",round(gc_total, 2)*100,"%"))
j <- 1
if (cu.bias){
cex = 0.7 * text.size
} else {
cex = 0.8 * text.size
}
for (i in which(c((length(ns) != 0), genome.length, total.gc))) {
graphics::text(0,0.05 - 0.05 * j, text[i],
font=4, cex = cex, col = info.color)
j <- j + 1
}
# Highlight IR
# Add legend
if (legend){
legend(x = -1.27, y = -0.6, legend = color_table$label, cex = 0.6 * text.size,
fill = color_table$col, bg = rgb(255, 255, 255, alpha = 0, maxColorValue = 255),
box.col = rgb(255, 255, 255, alpha = 0, maxColorValue = 255))
}
if (save){
grDevices::dev.off()
}
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.